直交表と応答曲面法の違い｜Pythonで実験計画と条件探索を実践

この記事で学べること

実験条件を少しずつ変えているのに、どの因子が効いているのか分からない。よい条件が見つかっても、その周辺をどう調べればよいか迷う。本記事は、このような材料研究者が実験計画を組み、結果を読み、次の確認実験へ進むための入門です。

• 因子・水準・応答、主効果・交互作用・曲率の意味を、条件表から理解できます。
• L9の読み方と、交互作用が別の列の効果と重なる「交絡」の注意点を学べます。
• L18、通常のL16、混合型L16を、水準数と調べたい関係から選べます。
• PythonでL18の条件表を作り、CSVを読み、水準別平均と主効果図を出力できます。
• 応答曲面を推定し、指定範囲内の候補条件を選び、確認実験につなげられます。

先に結論を示すと、直交表は条件の配置を定める設計、応答曲面法は実験・モデル化・条件探索を組み合わせる方法です。 両者は完全な対立関係ではありません。目的と推定する関係に合った点配置を選ぶことが重要です［1］［3］。

Pythonを初めて使う方は、まず第1〜4章で条件表と平均の意味を理解してください。その後、教材コードを実行し、出力CSVと図を確認すると読み進めやすくなります。実験計画法全体の位置づけは、材料研究者向けの実験計画法の基礎でも整理しています。

ここで使う薄膜熱処理の条件・応答は、すべて編集部の仮想例です。実測データ、公開研究の成果、推奨製造条件ではありません。 応答は単位のない教材指標で、大きいほどよいと仮定します。数値を実装置にそのまま設定せず、材料と装置の安全範囲を別途確認してください。

1．因子・水準・応答を、実験の言葉で理解する

入力条件と測定結果を分ける

薄膜の熱処理で、雰囲気、温度、保持時間、昇温速度を変えるとします。変更する項目が「因子」、各項目に設定する値や種類が「水準」、実験後に測定する品質指標が「応答」です。雰囲気A・Bのような種類の違いと、温度のような数値の違いを区別しましょう。

A・Bは雰囲気の識別名であり、大小関係を意味しません。コードでは種類として扱います。一方、温度の300と400の差には物理的な意味があります。ただし、最初の主効果比較では数値因子も水準ごとのカテゴリとして扱い、直線的な関係を初めから押し付けないようにします。

実務では、応答を一つ決めることも大切です。「品質がよい」だけでは解析できません。例えば透過率、膜厚のばらつき、抵抗率、欠陥率など、測定方法・単位・望ましい方向を先に決めます。複数指標がある場合は、一つの平均点にまとめる前に、最低限満たす制約と改善したい指標を分けて記録します。

主効果・交互作用・曲率の違い

主効果は、ほかの因子の条件を平均したときの、その因子の水準による違いです。温度が300℃の実験群と500℃の実験群の平均差が、その例です。適切な条件配置がなければ、温度の差に別の因子の偏りが混ざるため、実験表が重要になります。

交互作用は、ある因子の効果が別の因子の条件に依存することです。雰囲気Aでは高温で改善し、Bでは高温で悪化するなら、「高温がよい」という平均だけでは判断できません。次の数値は、意味を理解するための別の仮想例です。

温度ごとに平均すると、どちらも70です。しかし、雰囲気ごとの変化は逆方向です。主効果が小さいことと、その因子を変えても何も起こらないことは同じではありません。実務では、既知の反応機構や装置挙動から、強そうな交互作用を計画段階で挙げておきます。

曲率は、条件と応答の関係が直線から曲がることです。例えば低温60、中温80、高温65なら、中間付近によい条件がありそうです。ただし、3点の数値だけで真の最適温度は決まりません。測定誤差とほかの因子の影響を考慮し、追加実験とモデルで関係を調べます。

2．L9直交表は読み方の例。交互作用に注意する

9条件に何を配置しているのか

標準的なL9（3⁴）は、3水準の列を4列持ちます。Lは表の名称、9は基本の実験条件数です。以下はNISTに示された標準表と同じ水準配置です［1］。

数字は水準の番号です。第1列を温度に使えば、1・2・3を300・400・500℃に対応させます。どの列でも各水準が3回ずつ現れ、任意の2列では3×3の水準の組合せが1回ずつ現れます。これが、異なる因子の主効果を比較する際に役立つ配置です。

3因子3水準の全組合せは27条件ですが、L9の3列を使えば9条件になります。この削減には、全組合せにあった情報の一部を手放すという意味があります。「27回と同じことが9回で全部分かる」とは説明できません。

交互作用が特定の列と重なる

L9は表が小さく読みやすい一方、因子の割り付けに注意が必要な直交表です。 上の標準表では、第1・第2列の因子の交互作用は、第3・第4列の効果と交絡します。3水準どうしの交互作用には4自由度があり、その成分を第3・第4列の各2自由度と分離できません。この列関係は、表の配置に基づいて確認する必要があります。

例えば温度と時間を第1・第2列、昇温速度を第3列に入れると、温度×時間の影響の一部が昇温速度の効果に見える可能性があります。交絡とは、異なる原因を今回のデータから分けて推定できない状態です。これは測定を丁寧に行うだけでは解消しません［6］。

列の関係は「第1列×第2列の交互作用 → 第3列・第4列の効果と重なる」と覚えるだけで終えず、疑っている交互作用に必要な列を確保できるか確認します。3水準の交互作用を、一つの空き列だけで完全に扱えると考えないことも重要です。

4列をすべて因子で埋めると、主効果だけでも切片1個と各因子2自由度×4で、計9個のパラメータになります。9観測では残差自由度がありません。3列だけなら主効果モデルの残差自由度は2ですが、それをそのまま独立な実験誤差とみなすことはできません。省略した交互作用などが混ざるためです。

本記事では、L9を構造を理解する教材に使い、実践の主例は混合水準のL18にします。ただし、L9が常に使えないという意味ではありません。目的、想定する効果、割り付け、追加実験を含めて妥当性を判断する必要があります。

3．実務に使う直交表を選ぶ：L18・L16・L36

混合水準直交表とは

混合水準は、2水準・3水準・4水準など、異なる水準数の因子を同じ計画に含めることです。材料研究では、処理方式は2種類、温度は3段階、原料グレードは4種類というように、水準数がそろわない場面があります。そこで、因子構成に合う表を候補にします［1］。

上付きの数字は、その水準数の列がいくつあるかを示します。例えばL18の3⁷は、3水準の列が7列あるという意味です。これらは代表例で、L16の混合型には2⁹×4²などもあります［2］。「L16」とだけ書くと表を特定できないため、水準構成も記載してください。

通常のL16（2¹⁵）は二水準型です。交互作用が列と重なることもあるため、L9より条件数が多いという理由だけで自由な割り付けが可能になるわけではありません。二水準因子が少数で交互作用を調べたい場合は、完全実施要因計画も候補になります。

L18や混合型L16は、実際の因子の水準構成に合わせて検討する有用な選択肢です。ただし「材料実験は必ず混合直交表を使う」というルールではありません。禁止組合せが多い場合、推定したい特定のモデルがある場合、組成比の合計制約がある場合は、D最適計画や混合計画なども比較します。

L18の利点：混合水準の因子を、18条件で見通しよく比較できる

L18は、2水準と3水準の因子が混在し、まず各因子の主効果の傾向と追加調査の優先順位を知りたい場面で、有力な選択肢です。 標準表に水準を対応させることで、条件表の作成、結果の集計、研究者間での共有を進めやすくなります［1］。

• 実際の因子構成に合わせやすい：処理雰囲気はA・Bの2種類、温度・時間・昇温速度は各3段階という今回の課題を、一つの表で扱えます。
• 全組合せより少ない条件で主効果を調べられる：今回の4因子の全組合せは2×3×3×3＝54条件ですが、L18では18条件に配置できます。全交互作用まで同じように分かるわけではありませんが、初期調査の実験予算を抑え、確認実験や次段階の探索にも配分できます。
• ほかの因子の偏りを抑えて比較できる：任意の2列の水準の組合せが均等に現れます。主効果モデルの下では、ほかの割り付け因子の主効果と重ならずに、水準別平均を比較できる点が利点です。
• 3水準で中間の挙動も確認できる：低・中・高の平均を見ることで、単調な変化か、中間がよい傾向かを調べ、次に詳しく探索する範囲を考えられます。

実務では「L18で主効果の傾向を調べる→機構の知識と合わせて重要因子・疑わしい組合せを選ぶ→組合せの追加実験や応答曲面で詳しく調べる→独立に確認する」と進められます。この段階分けが、L18を活用する一つの方法です。

主効果の調査に活用し、重要な交互作用は次の計画で確かめる

L18の主効果比較が可能なことと、全交互作用を分離できることは別です。今回は雰囲気、温度、時間、昇温速度の4因子を使います。主効果モデルは切片を含め8パラメータ、18観測に対する残差自由度は10です。しかし、その残差には測定誤差だけでなく、省略した交互作用やモデルの不適合も入ります。

「空き列があるから、その列を誤差にすればよい」と自動判断しないでください。独立な誤差を調べるには、実験をやり直す反復が重要です。交互作用が研究上の中心課題なら、それを推定できる条件配置と反復を優先します［6］。

表を選ぶ前に、目的、因子の種類と水準数、強そうな交互作用、禁止条件、反復に使える予算を一枚に書きます。その上で、列の割り付けと推定可能なモデルを確認します。表の最大列数は、すべての列を使い切ることを推奨する値ではありません。

4．PythonでL18の計画・平均・主効果図を作る

教材の割り付けと実施順序

第1列に雰囲気、第2列に温度、第3列に時間、第4列に昇温速度を割り付けます。列を減らしても、残した列どうしの水準の組合せの均等性は保たれます。この教材では、列どうしの組合せ頻度をコードで確認してから条件表を出力します。

上表は標準の条件番号順です。実施順序とは別に管理します。同じ雰囲気を続けて実施すると、雰囲気と装置の時間変化が重なる可能性があります。コードは再現可能な乱数で実施順序を付けますが、雰囲気切り替えが難しい装置では、完全なランダム化が実行できない場合もあります。その場合は、分割区法など、制約に合った設計と解析を検討します。

環境を用意し、まず仮想データで実行する

Pythonの仮想環境を作れる方は、そこで次のライブラリを導入してください。JupyterやColabでも使えますが、下記はターミナルから実行する形です。完全版を orthogonal_rsm.py として保存し、同じフォルダで実行します。

python -m pip install numpy pandas matplotlib scipy statsmodels
python orthogonal_rsm.py

初回は実験を行わなくても動くように、条件表と仮想の測定結果を自動作成します。生成された l18_plan.csv は測定値が空欄の計画、l18_demo.csv は教材結果入りの表です。実験順序は run_order、条件の識別は run_id で保持します。教材結果を自分の実測値と混同しないよう、ファイル名を分けています。

水準別平均を計算して読む

短いコードでは、次のように結果を集計できます。温度の各水準に属する応答を平均し、その水準の観測数も一緒に出します。

import pandas as pd
results = pd.read_csv('doe_results/l18_demo.csv')
print(results.groupby('temp_C')['response'].agg(['mean', 'count']))

この割り付けでは、300℃の観測は6件あります。平均は「6件の応答を足して6で割る」計算です。例えば仮に60、62、64、66、68、70なら、平均65です。この6数値は計算説明用であり、配布データの実行結果とは別です。

主効果図では、横軸を水準、縦軸を平均にして点を線で結びます。温度の平均差が大きければ、今回の範囲で温度の影響を詳しく調べる候補になります。ただし、傾きの大きさは今回の水準幅に依存します。℃と分の変化量が違うため、線の傾きだけで物理的な感度を比較しません。

平均が最も高い水準を一つずつ選んで組み合わせても、その組合せが実測済みとは限りません。 また、交互作用があると期待した性能が出ないこともあります。候補の組合せを選んだら、独立に確認実験を行います。本教材の仮想応答にも交互作用を含めており、単純な主効果だけで完全に説明できるデータにはしていません。

自分のデータに置き換える

l18_plan.csv と後述の rsm_plan.csv をコピーし、実施した順序と実測応答を記入します。教材と同じ条件構成で解析する場合は、両方のファイルを指定してください。

python orthogonal_rsm.py --l18-input my_l18.csv --rsm-input my_rsm.csv --output my_results

入力の欠測、重複ID、非数値、想定外の条件配置は、ゼロで埋めず停止するようにしています。実験失敗は記録して原因を確認し、補完値で「測定したこと」にしません。計画どおりの入力であっても、測定法や材料の妥当性まで自動確認できるわけではありません。

この完全版は、教材の水準・条件数を検証する設計です。温度範囲や因子数を変える場合は、計画生成、対応表、入力検証、回帰モデルを一緒に変更する必要があります。CSVの数値だけを任意に変えて使える汎用ソフトではありません。記録項目の設計は、材料データの品質を高める方法も参考になります。

5．応答曲面法は、条件と応答の関係をモデルにする

一次項・交互作用・二次項を分ける

応答曲面法では、実験点を配置し、条件と応答の関係を回帰モデルなどで表し、よい条件や感度を調べます。二次モデルは曲率を表現する代表的な形ですが、すべての現象が二次式で十分に表せるわけではありません［3］。

ここでは雰囲気をB、昇温速度を5に固定し、温度と時間の2因子を扱います。これは説明を分かりやすくするために設定した第2段階です。L18の結果から、この固定条件の最適性が証明されたという意味ではありません。また、第2段階の教材応答は別に生成しており、L18との性能改善幅を比較する用途には使えません。

温度を x=(温度−400)/100、時間を y=(時間−20)/10 と変換すると、低・中・高が−1・0・1になります。これをコード化と呼び、数値の大きさをそろえ、係数と実験点を扱いやすくします。解析後は元の単位に戻します。

モデルは「応答の予測＝切片＋xの項＋yの項＋x×yの項＋x²の項＋y²の項」です。切片は中心条件の予測、一次項は中心付近の傾向、交互作用項は片方の効果がもう片方によって変わる関係、二次項は曲がりを表します。二次項があれば必ず範囲内に最大値があるとは限りません。

温度と時間を2水準で調べ、中心条件だけ追加しても、x²とy²の効果を別々に推定できません。中心点によって全体の曲率を検出できることと、因子ごとの二次項を推定できることを区別します［3］。

6．面中心型の中心複合計画を11回で組む

四隅・辺の中央・中心を配置する

本教材は、面中心型の中心複合計画を使います。2因子の四隅4条件、辺の中央4条件、中心条件3回で計11回です。異なる設定は9条件ですが、中心で3回の独立実験を行うため、観測数は11になります［4］。

面中心型は、星点が各因子の設定範囲内に収まるため、ここでは分かりやすい例として選びました。ただし回転可能な計画ではなく、中心から同じ距離でも方向によって予測の分散が変わります。別の中心複合計画では、星点が当初の上下限の外側に出るものもあり、設計の種類と安全範囲を確認します［4］。

中心3回は、同じ薄膜を3回測定するという意味ではありません。試料の準備や熱処理から独立に実施し、工程を含む変動を確認することを想定しています。同じ試料の繰り返し測定だけでは、試料作製や工程の変動を同じように評価できません。

自由度を確認し、モデルを詰め込みすぎない

2因子の二次モデルは切片を含め6パラメータです。この配置では設計行列の階数が6になり、11観測に対する残差自由度は5です。9つの異なる条件と中心反復を使い、純誤差の自由度2、不適合の自由度3に分けられます。コードでもこれらを確認します。

不適合とは、モデルの形が測定結果を十分に表せていない可能性です。中心の反復から得たばらつきと、モデルのずれを比較します。ただし、中心3回だけでは誤差推定の情報が少なく、等分散性や独立性が崩れると検定の解釈も変わります。有意にならなかったことは、モデルの正しさの証明ではありません。

7．Pythonで二次回帰・図・確認実験候補を出す

式を指定して回帰する

statsmodelsの式では、x:y は交互作用、I(x**2) は二次項です。式の意味を理解してから実行しましょう。公式資料にも、式を使ったモデルの指定方法が示されています［5］。

import statsmodels.formula.api as smf
rsm = pd.read_csv('doe_results/rsm_demo.csv')
model = smf.ols('response ~ x + y + x:y + I(x**2) + I(y**2)', data=rsm).fit()
print(model.params)

実測入力の場合は温度・時間からx・yを計算します。CSVに書かれたコード化値を無条件に信用せず、完全版では物理量から再計算します。係数を解釈するときは、コード化の基準、固定した雰囲気と昇温速度を併記してください。

残差と実施順序を見る

予測値と実測値の図では、点が対角線から大きく外れていないか見ます。残差は「実測−予測」で、予測値に対する曲がった傾向、ばらつきの広がり、特定の実験順での偏りを点検します。R²が高いだけで十分とは判断しません［6］。

本教材の残差図は学習に使った11件の結果です。新しい試料での予測性能を確認した図ではありません。日、原料ロット、測定器、洗浄や校正の記録と照合して、偏りの原因を調べます。異常点は、モデルに都合が悪いという理由だけで削除しないでください。

範囲内の候補を出し、確認実験を計画する

完全版では、x・yを−1〜1の各101点に分けた10,201点で予測し、平均予測が最大の格子点を候補にします。これは指定範囲と格子上での探索であり、連続空間の厳密な最適値や実材料の最大性能を保証しません。温度300〜500℃、時間10〜30分以外へ探索は広げません。

教材を実行すると、雰囲気B・昇温速度5℃/分の下で、候補は温度450℃、時間25分となりました。平均予測は87.210、平均応答の95％信頼区間は86.735〜87.684です。不適合検定のF値は6.030、p値は0.146でした。これは生成した仮想データの計算結果であり、p値だけでモデルの妥当性や条件の優位性を確定しません。

候補には平均応答の95％信頼区間も出力します。これは回帰モデルの仮定の下で、条件を固定した平均の不確かさを示すものです。次の一回の観測が入る範囲とは異なります。また、多数の候補から最大値を選んだ影響を調整した区間ではありません。区間が狭くても未考慮の機構変化や装置差は含まれません。

確認実験では、候補条件だけでなく、比較基準となる現行条件も同じ測定法で実施します。日やロットを変えた再現性、応答以外の品質制約、実施しやすさも確認します。候補が境界に来た場合は、範囲を広げる前に安全性とモデルの形を見直します。

8．使い分けと、よくある失敗

直交表で因子を絞って応答曲面へ進む流れは便利ですが、必須の順序ではありません。最初から2因子の曲率が主題なら、応答曲面計画から始める選択もあります。ベイズ最適化との目的の違いは、実験計画法とベイズ最適化の使い分けを参照してください。

失敗の一つは、平均だけを見て因子を外すことです。第1章のように、平均差が小さくても交互作用が強い場合があります。機構の知識と組合せごとの挙動を見て、除外する因子の根拠を残します。

もう一つは、計画表を順番どおりに実施し、装置変化と因子効果を重ねてしまうことです。ランダム化できない場合は、その制約を隠さず設計に組み込みます。反復、ブロック、分割区法の検討は、実施前に行うほうが有効です。

二次式への過信にも注意します。相変化、析出、分解、装置の切り替わりなどで関係が急に変わる場合は、一枚のなめらかな曲面が適切でないことがあります。範囲を分ける、追加因子を記録する、別のモデルを使うなど、現象と測定を見直します。

逐次探索を学びたい方は、Pythonで次の実験条件を選ぶベイズ最適化の記事へ進めます。本記事の格子探索は、既に得た回帰モデルから候補を選ぶ処理であり、予測不確かさを使って探索と活用を調整するベイズ最適化とは異なります。

9．実施前チェックリスト・FAQ

• 改善する応答、単位、測定法、望ましい方向を決めたか。
• 因子の水準数、実行可能範囲、禁止組合せを確認したか。
• 想定する交互作用と、列の交絡を確認したか。
• 表名だけでなく、L18（2¹×3⁷）などの水準構成を記録したか。
• モデルの階数、観測数、残差自由度を確認したか。
• 独立な反復と、同じ試料の再測定を区別したか。
• 実施順序、日、ロット、装置条件、失敗理由を残すか。
• 候補を独立実験で確認し、実装の判断責任者を決めたか。

L9は使わないほうがよいですか？

用途と割り付けが合えば候補になります。ただし小さい表ほど手放す情報を理解する必要があります。交互作用が別の列に現れる点を知らず、列を埋めて平均だけで結論づける使い方は避けます。

L18なら交互作用を気にしなくてよいですか？

L18は、混合水準の因子について主効果の傾向を整理するのに役立ちます。交互作用の確認は、目的に合わせて追加します。まず影響の大きそうな因子と組合せを選び、その交互作用を推定できる追加計画で調べましょう。特定の交互作用の推定が最初から主目的なら、完全実施要因計画なども候補にします。反復は誤差評価に役立ちますが、同じ配置の交絡の解消には条件配置の変更が必要です。

実験が失敗して欠測が出たら？

ゼロや全体平均を入れず、失敗原因と条件を記録します。再実施が可能か、欠測後も推定するモデルの階数が保たれるかを点検します。失敗しやすい条件を削除すると、探索範囲や適用範囲自体が変わる場合もあります。

配布コードを自分の材料にそのまま使えますか？

解析の流れは参考になりますが、教材は固定した水準・条件数を対象にしています。測定法、因子範囲、反復、制約、モデルは自分の課題に合わせて変更してください。コードが動くことと、実験計画が妥当であることは別に確認します。

純誤差や不適合の検定を必ず行うべきですか？

独立な反復と適切な仮定があれば有用ですが、検定だけでモデルを採用しません。反復数が少ない場合は検出力も限られます。残差、現象の知識、追加条件での確認を合わせて判断します。

まとめ：表の名前より「何を推定したいか」を先に決める

L9で条件表の読み方と交絡を学び、実践では水準構成に合ったL18やL16などを選ぶ。連続因子の曲率と局所的なよい条件を調べる段階では、応答曲面を推定できる配置を使う。この順に考えると、表を選んでから目的を合わせる迷いを減らせます。

Pythonは条件表、平均、モデル、図、候補の出力を助けます。一方、どの変動を再現するか、どの交互作用を調べるか、候補を実装してよいかは、材料と工程に基づいた判断が必要です。確認実験を計画の一部に入れ、推定と実測の違いを残して進めましょう。

完全版Pythonコード

実行確認：Python 3.12.14、NumPy 2.3.5、pandas 2.2.3、statsmodels 0.15.0。実行環境の版はsummary.jsonにも保存します。

実行すると計画CSV2本、仮想結果CSV2本、平均表、モデルの詳細、係数、確認候補、図3枚、実行環境と診断情報を出力します。教材は実験順序と独立反復を前提としていますが、コードから実際の独立性を確認することはできません。

"""Educational synthetic thin-film DOE; never use numeric settings as recipes."""
from pathlib import Path
import argparse, json, itertools, platform
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import statsmodels
import statsmodels.formula.api as smf
from scipy.stats import f as fdist

L9=np.array([[1,1,1,1],[1,2,2,2],[1,3,3,3],[2,1,2,3],[2,2,3,1],[2,3,1,2],[3,1,3,2],[3,2,1,3],[3,3,2,1]])
L18=np.array([
[1,1,1,1,1,1,1,1],[1,1,2,2,2,2,2,2],[1,1,3,3,3,3,3,3],
[1,2,1,1,2,2,3,3],[1,2,2,2,3,3,1,1],[1,2,3,3,1,1,2,2],
[1,3,1,2,1,3,2,3],[1,3,2,3,2,1,3,1],[1,3,3,1,3,2,1,2],
[2,1,1,3,3,2,2,1],[2,1,2,1,1,3,3,2],[2,1,3,2,2,1,1,3],
[2,2,1,2,3,1,3,2],[2,2,2,3,1,2,1,3],[2,2,3,1,2,3,2,1],
[2,3,1,3,2,3,1,2],[2,3,2,1,3,1,2,3],[2,3,3,2,1,2,3,1]])

def check_oa(a):
    for i,j in itertools.combinations(range(a.shape[1]),2):
        ui,uj=np.unique(a[:,i]),np.unique(a[:,j])
        counts=[np.sum((a[:,i]==u)&(a[:,j]==v)) for u in ui for v in uj]
        if len(set(counts))!=1 or min(counts)==0:
            raise ValueError('Pairwise balance failed')

def plans(out):
    check_oa(L9);check_oa(L18)
    rng=np.random.default_rng(42)
    d=pd.DataFrame({'run_id':[f'L18_{i:02d}' for i in range(1,19)],
        'atmosphere':np.where(L18[:,0]==1,'A','B'),
        'temp_C':np.array([300,400,500])[L18[:,1]-1],
        'time_min':np.array([10,20,30])[L18[:,2]-1],
        'ramp_C_min':np.array([2,5,8])[L18[:,3]-1]})
    d['run_order']=rng.permutation(np.arange(1,19))
    d['response']=np.nan
    d.to_csv(out/'l18_plan.csv',index=False)
    x=(d.temp_C-400)/100;y=(d.time_min-20)/10;z=(d.ramp_C_min-5)/3
    # Known synthetic generator for demonstration, not material measurements.
    demo=d.copy();demo['response']=70+4*(d.atmosphere=='B')+6*x-3*x*x+2*y-2*y*y+z+4*x*y+rng.normal(0,.8,18)
    demo.to_csv(out/'l18_demo.csv',index=False)
    pts=[(-1,-1),(-1,1),(1,-1),(1,1),(-1,0),(1,0),(0,-1),(0,1),(0,0),(0,0),(0,0)]
    r=pd.DataFrame(pts,columns=['x','y'])
    r['run_id']=[f'RSM_{i:02d}' for i in range(1,12)]
    r['atmosphere']='B';r['ramp_C_min']=5.
    r['temp_C']=400+100*r.x;r['time_min']=20+10*r.y
    r['run_order']=rng.permutation(np.arange(1,12));r['response']=np.nan
    r.to_csv(out/'rsm_plan.csv',index=False)
    rd=r.copy();rd['response']=85+5*r.x+3*r.y-6*r.x**2-4*r.y**2+2*r.x*r.y+rng.normal(0,.5,11)
    rd.to_csv(out/'rsm_demo.csv',index=False)
    return demo,rd

def validate(d,kind):
    required=['run_id','run_order','atmosphere','temp_C','time_min','ramp_C_min','response']
    if not set(required)<=set(d):raise ValueError('Missing columns')
    if d[required].isna().any().any():raise ValueError('Missing values; review failed runs, do not fill with zero')
    if d.run_id.duplicated().any() or d.run_order.duplicated().any():raise ValueError('Duplicate ID or run_order')
    for c in ['temp_C','time_min','ramp_C_min','response','run_order']:
        d[c]=pd.to_numeric(d[c],errors='raise')
        if not np.isfinite(d[c]).all():raise ValueError('Nonfinite numeric value')
    if kind=='l18':
        if len(d)!=18:raise ValueError('Expected all 18 planned observations')
        expected={'atmosphere':{'A','B'},'temp_C':{300,400,500},'time_min':{10,20,30},'ramp_C_min':{2,5,8}}
        for c,levels in expected.items():
            if set(d[c])!=levels:raise ValueError('Unexpected levels')
        a=d[['atmosphere','temp_C','time_min','ramp_C_min']].to_numpy()
        check_oa(a)
    else:
        if len(d)!=11 or set(d.atmosphere)!={'B'} or not np.allclose(d.ramp_C_min,5):raise ValueError('Expected fixed B and ramp 5, 11 runs')
        d['x']=(d.temp_C-400)/100;d['y']=(d.time_min-20)/10
        got=d.groupby(['x','y']).size().to_dict()
        exp={(a,b):1 for a,b in [(-1,-1),(-1,1),(1,-1),(1,1),(-1,0),(1,0),(0,-1),(0,1)]};exp[(0,0)]=3
        if got!=exp:raise ValueError('Expected face-centered design with 3 independent center runs')
    return d

def analyze_l18(d,out):
    d=validate(d.copy(),'l18');parts=[]
    fig,axes=plt.subplots(1,4,figsize=(12,3))
    for ax,c in zip(axes,['atmosphere','temp_C','time_min','ramp_C_min']):
        g=d.groupby(c).response.agg(['mean','count']);g['factor']=c;g['level']=g.index;parts.append(g.reset_index(drop=True))
        ax.plot(g.level.astype(str),g['mean'],'o-');ax.set_title(c);ax.set_ylabel('Synthetic response / a.u.')
    fig.tight_layout();fig.savefig(out/'l18_main_effects.png',dpi=130);plt.close(fig)
    pd.concat(parts).to_csv(out/'l18_level_means.csv',index=False)
    m=smf.ols('response ~ C(atmosphere)+C(temp_C)+C(time_min)+C(ramp_C_min)',data=d).fit()
    if np.linalg.matrix_rank(m.model.exog)!=8:raise ValueError('Main-effect model rank mismatch')
    (out/'l18_model.txt').write_text(m.summary().as_text())
    return {'n':len(d),'rank':8,'residual_df':int(m.df_resid),'note':'Residual includes omitted interactions; not pure error'}

def analyze_rsm(d,out):
    d=validate(d.copy(),'rsm')
    m=smf.ols('response ~ x+y+x:y+I(x**2)+I(y**2)',data=d).fit()
    rank=np.linalg.matrix_rank(m.model.exog)
    if rank!=6 or m.df_resid!=5:raise ValueError('Quadratic model rank mismatch')
    (out/'rsm_model.txt').write_text(m.summary().as_text())
    pd.DataFrame({'term':m.params.index,'coefficient':m.params.values}).to_csv(out/'rsm_coefficients.csv',index=False)
    pure=0.;dfpure=0
    for _,g in d.groupby(['x','y']):
        pure+=float(((g.response-g.response.mean())**2).sum());dfpure+=len(g)-1
    sslof=max(float(m.ssr)-pure,0.);dflof=int(m.df_resid-dfpure)
    F=(sslof/dflof)/(pure/dfpure) if pure>0 else None
    p=float(fdist.sf(F,dflof,dfpure)) if F is not None else None
    # Mean-response confidence interval, not future-observation interval.
    axis=np.linspace(-1,1,101);xx,yy=np.meshgrid(axis,axis)
    grid=pd.DataFrame({'x':xx.ravel(),'y':yy.ravel()})
    pred=m.get_prediction(grid).summary_frame(alpha=.05)
    ix=int(pred['mean'].to_numpy().argmax())
    candidate={'temp_C':float(400+100*grid.iloc[ix].x),'time_min':float(20+10*grid.iloc[ix].y),'atmosphere':'B','ramp_C_min':5.,'predicted_mean':float(pred.iloc[ix]['mean']),'mean_ci_low':float(pred.iloc[ix].mean_ci_lower),'mean_ci_high':float(pred.iloc[ix].mean_ci_upper)}
    pd.DataFrame([candidate]).to_csv(out/'confirmation_candidate.csv',index=False)
    fig,ax=plt.subplots(figsize=(6,4));cs=ax.contourf(400+100*xx,20+10*yy,pred['mean'].to_numpy().reshape(xx.shape),levels=15)
    ax.scatter(d.temp_C,d.time_min,c='black',s=20);ax.scatter(candidate['temp_C'],candidate['time_min'],marker='*',c='red',s=130)
    ax.set(xlabel='Temperature / C',ylabel='Time / min',title='Synthetic response: fixed atmosphere B, ramp 5');fig.colorbar(cs,ax=ax,label='Predicted mean / a.u.');fig.tight_layout();fig.savefig(out/'rsm_contour.png',dpi=130);plt.close(fig)
    fig,axes=plt.subplots(1,3,figsize=(11,3));axes[0].scatter(m.fittedvalues,d.response);lo=min(d.response.min(),m.fittedvalues.min());hi=max(d.response.max(),m.fittedvalues.max());axes[0].plot([lo,hi],[lo,hi],'k--');axes[0].set(xlabel='Fitted',ylabel='Observed')
    axes[1].scatter(m.fittedvalues,m.resid);axes[1].axhline(0,color='k');axes[1].set(xlabel='Fitted',ylabel='Residual')
    order=d.sort_values('run_order').index;axes[2].plot(d.loc[order,'run_order'],m.resid.loc[order],'o-');axes[2].axhline(0,color='k');axes[2].set(xlabel='Run order',ylabel='Residual')
    fig.tight_layout();fig.savefig(out/'rsm_diagnostics.png',dpi=130);plt.close(fig)
    return {'n':len(d),'rank':int(rank),'residual_df':5,'pure_error_df':dfpure,'lack_of_fit_df':dflof,'lack_of_fit_F':F,'lack_of_fit_p':p,'candidate':candidate}

def main():
    parser=argparse.ArgumentParser();parser.add_argument('--output',default='doe_results');parser.add_argument('--l18-input');parser.add_argument('--rsm-input');args=parser.parse_args()
    if bool(args.l18_input)!=bool(args.rsm_input):parser.error('Provide both input files, or neither for synthetic demo')
    out=Path(args.output);out.mkdir(parents=True,exist_ok=True)
    if args.l18_input:
        a=pd.read_csv(args.l18_input);b=pd.read_csv(args.rsm_input);mode='User supplied; physical validity not checked'
    else:a,b=plans(out);mode='Synthetic demonstration, not measured data'
    summary={'mode':mode,'versions':{'python':platform.python_version(),'numpy':np.__version__,'pandas':pd.__version__,'statsmodels':statsmodels.__version__},'l18':analyze_l18(a,out),'rsm':analyze_rsm(b,out)}
    (out/'summary.json').write_text(json.dumps(summary,ensure_ascii=False,indent=2));print(json.dumps(summary,ensure_ascii=False,indent=2))
if __name__=='__main__':main()

参考文献

［1］NIST/SEMATECH. e-Handbook of Statistical Methods, §5.3.3.10, Three-level, mixed-level and fractional factorial designs. 更新型Web資料、発表年不詳、2026年10月8日閲覧。https://www.itl.nist.gov/div898/handbook/pri/section3/pri33a.htm

［2］HBK ReliaSoft. Experiment Design and Analysis Reference: Taguchi Orthogonal Arrays. 更新型公式資料、発表年不詳、2026年10月8日閲覧。https://help.reliasoft.com/reference/experiment_design_and_analysis/doe/taguchi_orthogonal_arrays.html

［3］NIST/SEMATECH. e-Handbook of Statistical Methods, §5.3.3.6, Response surface designs. 更新型Web資料、発表年不詳、2026年10月8日閲覧。https://www.itl.nist.gov/div898/handbook/pri/section3/pri336.htm

［4］NIST/SEMATECH. e-Handbook of Statistical Methods, §5.3.3.6.1, Central Composite Designs. 更新型Web資料、発表年不詳、2026年10月8日閲覧。https://www.itl.nist.gov/div898/handbook/pri/section3/pri3361.htm

［5］statsmodels. Formulas: Fitting models using R-style formulas. 更新型公式ドキュメント、発表年不詳、2026年10月8日閲覧。https://www.statsmodels.org/stable/examples/notebooks/generated/formulas.html

［6］NIST/SEMATECH. e-Handbook of Statistical Methods, §5.4.3, How to model DOE data. 更新型Web資料、発表年不詳、2026年10月8日閲覧。https://www.itl.nist.gov/div898/handbook/pri/section4/pri43.htm

参考資料はいずれも公式Web資料で、個別ページのDOIは確認できないため記載していません。表の具体的な交絡、階数、自由度、教材の数値は、参考資料の考え方に加え、提示した配置を編集部が数値検証したものです。