Claude Codeで実験計画法の計画表をpyDOE3で作る手順

実験計画法の計画表は、因子と水準が決まれば機械的に作れます。Claude CodeにPythonパッケージのpyDOE3を使わせれば、フルファクトリアル、直交表、ラテン超方格の行列がその場で出ます。ただし「それらしい表」と「直交している表」は別物なので、生成のあとに交絡と直交性を検算するコードまで書かせるのが要点です。この記事では、手元で動かして確かめた出力を載せながら、計画表づくりの流れを示します。

pyDOE3はメンテナンスが終わっている

最初に押さえる点があります。pyDOE3のドキュメントのトップには、このパッケージは今後積極的にはメンテナンスされず、開発は本家のPyDOEに戻ったと書かれています。新機能、バグ修正、長期サポートはPyDOE側で提供されるため、PyDOEへの移行が勧められています。

それでもpyDOE3を使う場面はあります。既存の解析スクリプトがpyDOE3を前提にしている、あるいは本記事の出力と同じ結果を再現したい、といったときです。本記事の出力はPyPIにあるpyDOE3 1.6.2で取ったものです。新規に始めるなら、PyDOEのドキュメントで関数名と引数を確かめてからClaude Codeに渡すのが安全です。PyDOEの関数がpyDOE3と同じ名前・引数かどうかは、移行先のドキュメントで照らしてください。

どの関数で何の表が作れるか

用途ごとの対応は次のとおりです。ここを先にClaude Codeへ渡すと、関数選びの往復が減ります。

直交表の名前は list_orthogonal_arrays() で一覧できます。1.6.2では L4(2^3)、L8(2^7)、L9(3^4)、L12(2^11)、L16(2^15)、L18(6^1 3^6)、L27(2^1 3^12) など17種類が返りました。L9(3^4) は9回で3水準4因子の表です。

Claude Codeに渡す下準備

実行環境と約束事を先に固定します。仮想環境を作り、版を固定してインストールします。

python3 -m venv .venv
.venv/bin/pip install pyDOE3==1.6.2

CLAUDE.mdには、計画表づくりで守らせたいことを短く書きます。CLAUDE.mdは毎回の会話の冒頭で読み込まれる文脈であり、強制される設定ではありません。指示は具体的で簡潔なほど守られやすいと説明されています。

## 実験計画のルール
- 計画表は pyDOE3(.venv)で生成し、手計算で表を書かない
- 生成後に必ず検算する(交絡の一覧、全列ペアの水準組の数)
- 検算の出力をそのまま貼り、「直交している」と言うのは出力を読んでから
- 実験の実行順は乱数で並べ替えた列を別に持つ

毎回の実行確認を減らすなら、権限設定で該当コマンドだけを許可します。ワイルドカードは先頭のプログラム名とサブコマンドで範囲が決まるため、python 全般ではなく仮想環境のPythonに絞ります。

{
  "permissions": {
    "allow": ["Bash(.venv/bin/python *)"]
  }
}

ただし、これでは任意のスクリプトが走ります。計画表づくりの専用ディレクトリで使うか、doe.py 固定のルールに絞るかは、扱うデータの重さで決めます。

計画表を生成して検算するスクリプト

Claude Codeに「次の因子で計画表を作り、交絡と直交性を検算するスクリプトを書いて実行して」と頼むと、次のような形が得られます。以下は、実際に動かして確かめた形です。

import numpy as np
from itertools import combinations
from pyDOE3 import (fullfact, fracfact, fracfact_aliasing,
                    get_orthogonal_array, taguchi_design, lhs)
 
# 1. フルファクトリアル: 2水準 x 3水準 x 3水準
print("full:", fullfact([2, 3, 3]).shape)
 
# 2. 2水準4因子を8回に減らす(d = abc の交絡を許す)
d = fracfact("a b c abc")
print("runs:", d.shape)
for line in fracfact_aliasing(d)[0]:
    print(line)
 
# 3. 直交表L9の検算: 全列ペアで水準の組が同数ずつ出るか
oa = get_orthogonal_array("L9(3^4)")
for i, j in combinations(range(oa.shape[1]), 2):
    _, counts = np.unique(oa[:, [i, j]], axis=0, return_counts=True)
    print(i, j, "組の数:", len(counts), "各組の回数:", set(counts.tolist()))

出力は次のとおりです(pyDOE3 1.6.2)。

full: (18, 3)
runs: (8, 4)
a = bcd
b = acd
c = abd
d = abc
ab = cd
ac = bd
ad = bc
abcd
0 1 組の数: 9 各組の回数: {1}
0 2 組の数: 9 各組の回数: {1}
0 3 組の数: 9 各組の回数: {1}
1 2 組の数: 9 各組の回数: {1}
1 3 組の数: 9 各組の回数: {1}
2 3 組の数: 9 各組の回数: {1}

読み方を順に示します。

• フルファクトリアルは2×3×3で18行。回数が多すぎるなら次の一部実施へ進みます
• fracfact("a b c abc") は4因子を8回で回す式で、4列目を abc の積で作ります
• 交絡の一覧では、主効果 a が3因子交互作用 bcd と区別できません。一方で2因子交互作用は ab = cd のように対になって混ざります
• 最小の別名が4項(abcd)なので、ドキュメントの定義に従えば分解能IVです。主効果は2因子交互作用と混ざりません
• L9は6組の列ペアすべてで、水準の9通りの組が1回ずつ現れました。これが直交の確認です

主効果だけを見るなら、この8回で足ります。2因子交互作用を個別に見たいときは、ab と cd を分けられないため、折り返し(fold)で追加の8回を足すか、回数の多い表を選びます。

水準を実際の値に置き換える

直交表の中身は0始まりの水準番号です。taguchi_design に実際の水準を渡すと、設定値の表になります。

levels = [["low", "mid", "high"], [100, 200, 300],
          ["A", "B", "C"], [1, 2, 3]]
print(taguchi_design("L9(3^4)", levels)[:3])

[['low' 100 'A' 1]
 ['low' 200 'B' 2]
 ['low' 300 'C' 3]]

表の列の順は因子の順です。置換後の表は、実行順の列を足してから実験担当者に渡します。実行順は表の並びのまま回すと、時間経過による変動が特定の因子に乗るため、乱数で並べ替えた列を持たせます。

連続値の因子はラテン超方格にする

温度や圧力のように連続的に動かせる因子で、水準を決め打ちしたくないときは lhs を使います。ここでの落とし穴は、ラテン超方格は直交表と違い、列同士の相関が0になる保証がないことです。

手元で lhs(3, samples=8, criterion="maximin", seed=42) を実行し、列間の相関係数を取ると、絶対値の最大が0.43でした。criterion="correlation" に変えても0.41で、この条件ではほとんど下がっていません。8点・3因子では、criterionを替えれば相関が消えるとは限らないので、相関行列を出して見る検算を必ず入れます。seed 引数は1.6.2で通りましたが、ドキュメントの lhs の書式には載っていません。

点数の既定値は因子数と同じです。実験回数に合わせて samples を必ず指定します。分布を変えたいときは、0から1の値を scipy.stats の逆関数(ppf)で目的の分布に変換します。

x = lhs(3, samples=8, criterion="maximin", seed=42)
r = np.corrcoef(x.T)
print(round(float(np.abs(r - np.eye(3)).max()), 2))

ドキュメントの例をそのまま写さない

pyDOE3のドキュメントの例には、Python 3では動かないものが混じっています。lhs の節の分布変換の例は xrange を使っており、Python 3では range に直さないとエラーになります。random_k_means の例では random_state が使われていますが、同じ節の引数一覧では非推奨で、seed を使うよう書かれています。

ほかに、コードを書かせる前に伝えておくと手戻りが減る仕様が3つあります。

• fracfact の生成式は小文字か大文字のどちらかに統一します。混在はできません
• fullfact は水準番号を0始まりで返し、ff2n と fracfact は-1と1で返します。コードの意味づけを揃えないと、解析で水準を取り違えます
• get_orthogonal_array の水準も0始まりです

Claude Codeには、ドキュメントのコード例をそのままコピーさせず、手元で一度実行して動いた形だけを残させます。検算の出力が空でないことも見ます。

回数を減らす・表を選ぶときの判断

水準数と因子数が決まったら、次の順で表を選びます。

pbdesign は4の倍数の回数で作る2水準の表です。pbdesign(7) は8行、pbdesign(8) は12行、pbdesign(11) も12行で返りました。7因子8回の pbdesign(7) と fracfact("a b ab c ac bc abc") は、行を並べ替えれば同じ集合になります。ドキュメントは両者が同一だと説明し、np.all(pbdesign(7)==fracfact(...)) が True になる例を載せていますが、1.6.2では行の並びが違うため、この比較は False になりました。行の集合で比べるか、検算を全列ペアの数え上げに任せてください。

冗長な計算を避けたい場合は、生成式を探索させる fracfact_opt もあります。fracfact_opt(6, 2) は生成式として a b c d bcd acd を返しました。

タグチ流のSN比は compute_snr で計算できます。[90, 95, 93] を望大特性で渡すと、1.6.2は39.33 dBを返しました。ドキュメントの例には39.133...と書かれていますが、-10*log10(mean(1/y^2)) を手で計算すると39.33になるため、手元の出力と一致します。数値を報告書に載せるときは、ドキュメントの例ではなく手元の実行結果を使います。

交絡が大きすぎたときの追加実験

分解能IIIの表では、主効果が2因子交互作用と混ざります。fracfact("a b ab") は3因子を4回で回す表で、別名の一覧は a = bc、b = ac、c = ab でした。主効果を読みたいのに、その因子が他の因子の交互作用と区別できない状態です。

この場合は fold で全列の水準を反転した4行を足します。fold(fracfact("a b ab")) は8行を返し、ドキュメントでは、全列を反転すると主効果と2因子交互作用の連鎖が切れると説明されています。特定の列だけ反転するなら columns=[2] のように列番号を0始まりで渡します。回数は倍になるので、最初から分解能IVの表を選ぶ場合と比べて、追加実験で済ませるメリットがあるかをClaude Codeに表で比べさせます。

Claude Codeへの頼み方

一度に全部を頼むより、因子の表を渡して段階を分けるほうが検算が効きます。たとえば次のように進めます。

因子は温度(3水準)、時間(3水準)、圧力(3水準)、触媒(3種類)です。
1. pyDOE3のL9(3^4)で計画表を作り、実際の水準に置き換えてください
2. 全列ペアの水準の組が1回ずつ現れるか数え、結果を貼ってください
3. 実行順を乱数で並べ替えた列を足し、CSVに保存してください
4. CSVの行数と列数も確認してください

2の出力が空だったり、組の数が9未満だったりしたら、表の取り違えか列の読み違いです。割付の説明を文章で返させるだけの依頼は、検算にならないので避けます。

次の工程

計画表ができたら、実験結果から分散分析まで進めます。L8への割付、誤差列の見方、分散分析の手順はClaudeで実験計画法の直交表を組み、分散分析まで進める手順にあります。結果の差が偶然かどうかを2群で見るだけなら、ClaudeにA/Bテストの有意差を検定させる手順のほうが軽く済みます。

検算を毎回忘れずに走らせたい場合は、完了条件をテストの形で持たせる方法もあります。Claude Codeのagent hookで「テストが通るまで止めない」を検証させるの考え方を、直交性の検算スクリプトに当てはめられます。

まとめ

pyDOE3は関数1つで計画表を作れますが、メンテナンスは終わっており、新規ならPyDOEの確認が先です。表の出来は、交絡の一覧と全列ペアの水準組の数という2つの検算で判断します。連続値のラテン超方格は相関が残るため、直交表と同じ安心は得られません。Claude Codeには、生成と検算を同じスクリプトに入れさせ、出力を読ませてから結論を言わせます。