MATLAB で実験計画法（DOE​）と回帰分析を用いて​、総当たり実験よりも​コンパクトに説明変数​と応答変数を関連付け​るにはどうすればよい​ですか。

複数の因子が 1 つまたは複数の応答変数にどのような影響を与えるかを評価したいと考えていますが、各因子に複数の水準がある場合、総当たり実験では実験数が大きくなりすぎます。たとえば、4 つの因子をそれぞれ 3 水準で評価する場合、候補となる組合せ数はすぐに増えてしまいます。そこで、どの因子が重要かを見積もりつつ、よりコンパクトなワークフローを取りたいと考えています。

また、MATLAB 上で実験計画の段階と解析の段階をつなげたいと考えています。応答が単一系列の場合は線形回帰で因子の寄与度を推定したい一方、応答が複数系列の場合は、複数出力を同時に扱える多変量手法を使いたいと考えています。

Answer 1

[MathWorks Support Team]

実験回数を削減するために田口の直交表を用い、その後、応答データの形に応じて重回帰分析または部分最小二乗法（PLS）を適用します。

基本的な考え方は、より大きな総当たり実験の代わりに、因子列どうしの独立性を保つ直交計画を用いることです。これにより、より少ない実験回数でデータを取得しながら、各因子の相対的な影響度を見積もることができます。

たとえば、L9 直交表は 3 水準の 4 因子までを 9 回の実験で表現できます。MATLAB では、次のように実験計画を作成できます。

>> taguchiTypes
TypeSupportedFactors 
________________________________________________
L43 two-level factors 
L8 7 two-level factors 
4 two-level factors and 1 four-level factor 
L9 4 three-level factors 
L12     11 two-level factors 
…
…

>> dtag = taguchiDOE([-1 0 1],[-1 0 1],[-1 0 1],[-1 0 1],'L9');
>> dtag.Design
ans =
9×4 table
Factor1 Factor2 Factor3 Factor4
_______ _______ _______ _______
-1-1 -1 -1 
-1  0  0  0 
-1  1  1  1 
 0 -1  0  1 
 0  0  1 -1 
 0  1 -1  0 
 1 -1  1  0 
 1  0 -1  1 
 1  1  0 -1 

この種の直交表では、各因子列どうしの内積がゼロになるように構成されています。つまり、因子が高い独立性を持つ組合せでサンプリングされるため、より詳細な物理モデルや機構解析に進む前のスクリーニングや寄与度分析に適しています。

Example: Factor2 × Factor3
(-1)×(-1) + 0×0 + 1×1 + (-1)×0 + 0×(-1)
+ 1×(-1) + (-1)×1 + 0 ×(-1) + 1×0 
= 0

実際の例として、3 因子をそれぞれ 3 水準で評価する場合、L9 直交表を用いて 9 回の実験を定義できます。

>> dtag = taguchiDOE([200 400 600], …
[0 1 2],…
[0.1 0.5 1.0], 'L9');
>> T = dtag.Design
T =
  9×3 table
    Factor1    Factor2    Factor3
    _______    _______    _______

      200         0         0.1  
      200         1         0.5  
      200         2           1  
      400         0         0.5  
      400         1           1  
      400         2         0.1  
      600         0           1  
      600         1         0.1  
      600         2         0.5 

1. 実験結果が単一応答系列(最小二乗法を用いた多重回帰分析)の場合

実験結果が単一の応答ベクトルである場合は、最小二乗法による重回帰分析で因子の寄与度を推定できます。まずは、次のように直接フィットできます。

e.g.
>> a = T \ Y
a =
    0.0001
    0.1283
    0.3693

ただし、因子ごとにスケールが異なる可能性があるため、係数の大きさを比較する前に設計行列を正規化する必要があります。

e.g.
>> a = normalize(T) \ Y
a =
    0.1346
    0.0836
    0.0104

* Y : zero as the average, 1 as the standard deviation

正規化後、係数が大きいほど影響が大きいことを示します。これは、応答列が一つで単純な線形寄与モデルを求める場合に適したワークフローです。

2. 実験結果に2つ以上の応答系列が含まれている場合(部分最小二乗法)

正規化後は、係数が大きいほど、その因子の影響が大きいことを示します。これは、応答が単一系列で、シンプルな線形寄与モデルを使いたい場合に適したワークフローです。

実験結果が複数の応答系列を含む場合は、PLS に切り替えることができます。

e.g.
>> ncomp = 3; % Number of Principal Components
>> [XL, YL, XS, YS, beta, PCTVAR] = plsregress(T, Y, ncomp);
>> beta 
beta =
    0.1052
    0.0010
    0.0791
    0.0044

このモデルでは、`beta(1)` が切片、`beta(2:end)` が各因子の寄与度に対応します。予測値と残差は次のように計算できます。

e.g.
>> T_aug = [ones(size(T,1),1) T]; % intercept
>> Yhat = T_aug * beta;           % response（9×2）
Yhat =
    0.2786
    0.3150
    0.3920
    …

>> residuals = Y – Yhat;
>> stem(residuals)

これにより、実測値とモデルによる予測値を比較し、残差の挙動を確認できます。

PLS は応答側が多変量である場合に特に有効であり、非常に高次元な問題にも適用できます。一方で、オーバーフィットを避けるためには、潜在成分数を適切に選ぶ必要があります。実務上は、累積寄与率を確認する方法が有効です。

>> ncomp = 10; % Number of Principal Components
>> [XL, YL, XS, YS, beta, PCTVAR] = plsregress(T, Y, ncomp);
>> plot(1:10,cumsum(100*PCTVAR(2,:)),'-bo');

要約すると、ワークフローは次のとおりです。

1. L9 や L18 などの直交実験表を作成する。

2. 削減した実験セットを実施して応答を取得する。

3. 応答が単一の場合は、正規化した重回帰分析を用いる。

4. 応答が複数の場合は、PLS を用いる。

5. 因子の重要度を解釈する前に、係数、予測値、残差、および主成分数を確認する。

この方法により、総当たり実験を行わなくても、MATLAB 上で説明変数と応答変数を関連付けるための、コンパクトな実験・解析フローを構築できます。