PythonによるD最適化基準を用いた実験計画

はじめに

製造業に勤めている皆様の中には、実験計画法について学んだことが多いのではないでしょうか？（私もその中の一人でした。）
おそらく皆様が学んだ実験計画法では直交表に基づき、実験計画を立てる方法が記載されていたかと思います。
一方で、目標とする物性を達成するための条件探索という目的に従うのであれば、直交表対比で効率が良いD最適化基準による実験計画があります。
本記事ではPythonを使用したD最適化基準の実装を紹介します。
また、実験計画は製造業だけでなく、例えばマーケティング的な実験を行う際の水準選択としても有用かと思いますので、その関係の方にも本記事を参考にしていただければと思います。

D最適化基準とは？

実験計画法の目的

D最適化基準を説明する前に、実験計画法を使う目的について説明します。
実験計画法を使う目的は、なるべく少ない実験回数でなるべく多くの情報を得ることにあります。
これを達成するためには、各条件が独立している、つまり、似たような条件を繰りかえさないことが必要です。
似たような実験を繰り返さないこと、というのは極端な例としてある条件を等倍した実験を行うこと（ex. 温度、圧力ともに2倍してしまったらどちらが物性に効いたかわからない）が挙げられます。
我々が良く知る？直行表に基づき、実験を行えば何も考えなくとも上述した内容が考慮されて水準を組むことができます。
一方で、直交表では固定された数の要因しか考慮にいれることができず、柔軟性に欠けます。
そこで、直交表に代わる水準の決め方として、D最適化基準があります。

D最適化基準の概要

ある特徴量行列$X$が存在するとします。
尚、特徴量行列とは例えば以下のような実験水準を想定してください。

温度 (°C)	ポリマーA含有量 (g)	圧力 (kPa)
25	10	101.3
30	15	105.2
35	20	108.9
40	12	102.5
45	18	107.1

D最適化基準とは、端的に述べると、特徴量行列$X$とその転置行列である$X'$があるとき、以下式で与えられます。

$$Max(det(X'X))$$

つまり、特徴量行列とその転置行列を掛け合わせたものの行列式が最大化すれば良いということです。
この式が何を意味するのか考えていきましょう。
まず、$X'X$に関しては、各特徴量同士の内積が各成分に出現します。各特徴量同士が相関している場合、$X'X$の逆行列$(X^{\prime} X)^{-1}$は計算できない、つまり、行列式$det(X'X)$は小さくなります。
よって、逆に考えると、行列式$det(X'X)$が大きいほど各特徴量同士は相関していない、つまり情報量が最大化していると考えることができます。
また、別の捉え方として、$X'X$は線形回帰モデルの回帰係数を算出する際に出現します。

$$\hat{\beta} = (X^{\prime} X)^{-1} X^{\prime} y$$

上式より、逆行列$(X^{\prime} X)^{-1}$が計算できない、つまり、行列式$det(X'X)$が小さいと回帰係数は算出できません。
以上より、D最適化基準についてイメージを持って頂ければ幸いです。

PythonによるD最適化基準の実装

ここからは具体的にD最適化基準の計算に移りたいと思います。
まずは、必要な関数を定義します。
尚、上述していない点として、計算を安定させるために各列を正規化する処理を挟んでいます。

def normalized(X):
    # 各列の平均と標準偏差を計算
    mean = np.mean(X, axis=0)
    std_dev = np.std(X, axis=0)

    # 非常に小さい標準偏差で除算を避けるための小さい定数（例えば1e-8）を追加
    std_dev += 1e-8

    # 正規化
    X_normalized = (X - mean) / std_dev
    
    return X_normalized

def calculate_d(X_normalized):
    # X の転置を計算
    X_transpose = X_normalized.T

    # 情報行列 X'X を計算
    information_matrix = np.dot(X_transpose, X_normalized)
    
    # 行列式を計算
    d = np.linalg.det(information_matrix)
    
    return d

def main(X):
    X_normalized = normalized(X)
    d = calculate_d(X_normalized)
    return d

続いて、各特徴量が相関している行列（多重共線性を持つ行列）とそうでない行列を作成します。

# 多重共線性がある行列の作成
# ここでは、第2列は第1列のほぼ同じ値で、第3列は第1列と第2列の平均として作成されます。
X_multicollinearity = np.array([
    [2, 2.1, 2.05],
    [4, 4.2, 4.1],
    [6, 5.9, 6.0],
    [8, 8.1, 8.05]
])

# 多重共線性がない行列の作成
# ここでは、各列はランダムな値を持っており、列間の相関は低いです。
np.random.seed(0)  # 一貫した結果を得るためのシード値
X_no_multicollinearity = np.random.rand(4, 3)

これらの行列に対してmain関数を適用するとD最適化基準値が出力されます。
X_multicollinearityのD最適化基準:2e-06
X_no_multicollinearityのD最適化基準:26
出力結果より、列間の相関が高いとD最適化基準値がほぼ0になることがわかります。

実際に水準を決める際には、各列の上下限等の制約条件をかけた上で、水準をランダムに作成し、D最適化基準値が高い水準を選択すれば良いでしょう。

おわりに

本記事ではD最適化基準について紹介しました。
理論はやや複雑なところがあるかもしれませんが、Pythonによる実装は簡単であることがわかるかと思います。
一方で、実務上では既に実験の積み上げがある、強力な理論がある、実験コストがかかりすぎる、といったことから、上記で書いたような網羅的にパラメータを探索するような実験をする機会は少ないかもしれません。
ただ、開発初期・事業初期段階でパラメータの効き具合が掴めていないときは、本手法のような実験計画に基づき、水準を組むことは、後の人に知見を残すという意味でも行った方がいいかと思います。