中心複合計画と応答曲面法

中心複合計画は応答曲面法のための実験です。今回は実験結果から応答曲面法のパラメータ計算を追ってみます。

応答曲面法

応答曲面（Response Surface）とは、予測変数（Predictor variables）から応答（Response）ｙを関係近似したものです。[math] \displaystyle \varepsilon [/math]を誤差とすると下記のように表せます。

[math] \displaystyle y = f( x_1 \cdots x_n ) + \varepsilon [/math]

応答曲面法において関数の形に制限はないのですが、取り扱いが簡単な多項式近似が用いられることが多いです。多項式近似の他には、クリギング法、ニューラルネットワークやその派生のRBFネットワーク、スプラインなどの補間法も考えることができます。ただし左記のような非線形手法を適用すると、P値など統計的な評価が困難という問題があります。そして大抵の場合、少なくとも最適化対象範囲内では多項式近似で問題ないことが多いです。そんなわけで多項式近似を見ていきます。

最小二乗法

多項式近似のパラメータは最小二乗法で簡単に求めることができます。モデルは下記で、以下データの数を[math] \displaystyle n [/math]、近似パラメータの数を[math] \displaystyle k [/math]とします。

[math] \displaystyle y ={\beta}_0 + \sum_{i=1}^{n}{{\beta}_i x_i} + \sum_{i<j}^{n}{{\beta}_{ij} {x_i}{x_j}} + \sum_{i=1}^{n}{{\beta}_{ii} {x_i}^2} [/math]

2変数の場合は下記の通り。

[math] \displaystyle y ={\beta}_0 + {\beta}_1 x_1 + {\beta}_2 x_2 + {\beta}_{3} {x_1}{x_2} + {\beta}_{4} {x_1}^2 + {\beta}_{5} {x_2}^2 [/math]

これを行列式で表します。

[math] \mathbf{y}=\left \{ \begin{matrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{matrix} \right \} [/math]

[math] \mathbf{X}=\left \{ \begin{matrix} 1 & x_{11} & x_{12} & \cdots & x_{1n} \\ 1 & x_{21} & x_{22} & \cdots & x_{2n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & x_{n2} & \cdots & x_{nn} \\ \end{matrix} \right \} [/math]

[math] \mathbf{\beta}=\left \{ \begin{matrix} \beta_1 \\ \beta_2 \\ \vdots \\ \beta_n \end{matrix} \right \} [/math]

[math] \mathbf{\varepsilon}=\left \{ \begin{matrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{matrix} \right \} [/math]

と置くと、近似モデルの行列式表示は下記のようになります。

[math] \displaystyle \mathbf{y}= \mathbf{X} \beta + \varepsilon [/math]

最小二乗法から[math] \displaystyle \mathbf{\beta} [/math]の推定値[math] \displaystyle \mathbf{b} [/math]が下記のように得られます。

[math] \displaystyle \mathbf{b}= (\mathbf{X}^{\mathit{T}} \mathbf{X})^{-1} \mathbf{X}^{\mathit{T}} \mathbf{y} [/math]

残差平方和 [math] \displaystyle S_E [/math] (Square Sum of Errors)

[math] \displaystyle S_E= \mathbf{y}^{\mathit{T}} \mathbf{y} + \mathbf{b}^{\mathit{T}} \mathbf{X}^{\mathit{T}} \mathbf{y} [/math]

回帰二乗和 [math] \displaystyle S_R [/math] (Regression Sum of Squares)

[math] \displaystyle S_R = \mathbf{b}^{\mathit{T}} \mathbf{X}^{\mathit{T}} \mathbf{y} - \frac{1}{n} (\sum_{i=1}^{n}{y_i})^2 [/math]

全変動=応答ｙの平均値とデータの差の平方和 [math] \displaystyle S_T [/math] (total sum of squares)

[math] \displaystyle S_T = \mathbf{y}^{\mathit{T}} \mathbf{y} - \frac{1}{n} (\sum_{i=1}^{n}{y_i})^2 [/math]

決定係数[math] \displaystyle R^2 [/math]は

[math] \displaystyle R^2 = \frac{S_R}{S_T} =1 - \frac{S_E}{S_T} [/math]

決定係数については、説明変数/次数を多くとれば決定係数は大きくなります。このため、単位自由度あたりの残差を比較すると公正な比較になり、一般には自由度調整済み決定係数 [math] \displaystyle {R^2}_{adj} [/math]を計算します。

[math] \displaystyle {R^2}_{adj} =1 - \frac{S_E/(n-k-1)}{S_T/(n-1)} [/math]

AIC (Akaike's Information Criterion)

多項式近似を行うとパラメータ数が多いほど適合度(フィッティング精度)は上がりますが、過適合(overfitting)になる可能性があります。パラメータ数と適合度のバランスをとる指標にAIC (Akaike's Information Criterion)があります。具体的にはAICを最小化すれば良いモデルが得られることになります。応答曲面法の解析でもAICを計算しますが、式は下記のようになります。

[math] \displaystyle AIC = n{\mathrm{ln}} \frac{S_E}{n} + 2k + n{\mathrm{ln}} 2\pi [/math]

いくつか流儀があって、右辺第三項は定数項なので省いて下記で計算する場合もあります。

[math] \displaystyle AIC = n{\mathrm{ln}} \frac{S_E}{n} + 2k [/math]