応答曲面の観点からみた最適計画の比較
M2012MM031 奥村和也
指導教員：松田眞一

1

はじめに

3.1

基本概念

最適計画はモデルを基に計画を評価する基準を定めて， X を計画行列とすると，正規方程式で求められる応答その基準で最適化する計画である．1959 年に最適計画は 曲面の各項の係数の分散は，応答 y の分散 σ と正方行列
T
初めて提案されたが，多くの繰り返し計算を必要とし，計 X X の逆行列からなる．X ，y は，実験点数を n とする算が困難であることから，当初は理論的な面での研究しか と，次のように定義される．

あまりされなかった．それでも，1990 年半ば頃より，計1 x11 x12 · · · x1n
算機器の発達にともない，最適計画は実用化されるよう
 1 x21 x22 · · · x2n になってきた．そして現在では，計画の良し悪しを評価
(2)
X=
..
..
.. ..
 ..
.
.
.
. 
 .
する基準として，A 最適基準や D 最適基準など多くの最1 xn1 xn2 · · · xnn
適基準が用いられている．しかし，何を目的にするかで，計画を最適化する基準が変わることから，どの最適基準
が一番よいかということは一概には言えない．そこで今y1回の研究では，応答曲面の観点からどの計画が最適であ y2 るかについての比較を行う．
(3)y= . 
 .. 2 応答曲面について
yn
応答曲面法とは，予測される応答 y と予測変数 x1 , …
この分散を小さくするには σ を小さくすることも重要, xn の関係について，実験計画にしたがってデータを収T−1集し，それを解析することで応答と量的な因子の関係を であるが，(X X) の対角成分を大きくすれば係数の分散は相対的に小さくすることができる．この (X T X)−1探索する一連の方法である．
は近似モデルと実験点の座標だけで決定される．X T X を
y = a1 x1 + a2 x2 + … + an xn
実験点数 n で割ったものは情報行列と呼ばれている．

+a11 x21 + a22 x22 + … + ann x2n
M (ξ) =+a12 x1 x2 + a13 x1 x3 + … + an−1n xn−1 xn + ε(1)ここで，ε は誤差と呼ばれる．ai xi は因子 xi の１次の効果を，aii x2i は因子 xi の２次の効果を表し，aij xi xj は因子 xi と xj の１次 × １次の交互作用を表している．応答曲面法において関数の形に特に制限はないが，一般的には１次項，２次項，積項までである．応答曲面法は，製品プロセスの最適化やばらつきの減少などの品質工学の分野において，特にアメリカで実用化されている手法である．(山田 [1]，轟 [2] 参照)

3

計算機支援計画

電子計算機の発展に伴い，表を用いた実験計画から計算機支援の実験計画が注目されるようになってきた．計算機支援の実験計画は，あらかじめ実験計画の候補点をいくつか用意しておき，その中から次に述べる最適性を用いて最適実験計画の組み合わせを求めるものである．計
算機支援の実験計画は計算機の発展によって最近用いられるようになった．計算機支援の方法は使用法が簡単であり，現在は多く用いられている．計算機支援実験計画
はいくつか提案されており，アルファベットを用いて命名されている．A，D，I 最適基準などが有名であり，下記に説明する．

XT X
n(4)

ξ は連続計画であり，n 個の異なる点での試行を持つな
らば，wi を結合させた重みとして，次のように表すことができる．
[]x1 x2 · · · xn
ξ=(5)w1 w2 · · · wn
(山田 [1] 参照)

4

最適基準の種類

4.1

A 最適基準A 最適基準は応答曲面の係数の分散を最小化する最適化基準であり，情報行列の逆行列の対角成分の和 (trace)に着目し，これを最小化する実験計画の最適基準である．
min[tr(M −1 (ξ))](6)

A 最適基準は対角成分だけを考慮し，応答曲面の各項の係数の共分散は考慮しない．このため，A 最適な実験計画でも共分散の絶対値が大きい場合が存在し得る．共分散の絶対値が大きいことはその変数 xi と xj の実験計画に偏りがあることを意味する．つまり，回転可能性という観点からは良好ではない実験計画になり得る．(山田[1] 参照)

4.2

D 最適基準D 最適基準は情報行列の行列式を最大化する最適化基準である．
max |M (ξ)|(7)

変数の座標範囲が −1〜1 に正規化されている場合，D最適性の優劣を表す Def f (Deﬃciency) は次の式で表すことができる．

(|X T X|)1/p
Def f =n(8)ここで，p は未知係数の数である．Def f は最大値が１となり，大きいほど良い実験計画となる．
一般に逆行列の各成分にはその行列の行列式の逆数が掛けられる．D 最適基準は (X T X) 行列の行列式を最大化することで (X T X)−1 行列の全成分を相対的に最小化している．これによって，D 最適基準は情報行列の対角成分だけでなく，共分散成分も含んでいる．つまり，D最適基準は係数全体の分散を考慮した一般化分散を最小化させている．このため，D 最適基準で選択された点は回転可能性という点からも良好な実験計画である．
D 最適計画は固有である必要はない．もし ξ1∗ と ξ2∗ が最適計画であるならば，計画ξ ∗ = cξ1∗ + (1 − c)ξ2∗(0 ≤ c ≤ 1)

(9)もまた D 最適である．(Atkinson and Donev[3] 参照)

4.3

5.2 実行結果
計画の水準をそれぞれ５，６，４とする．
> levels.design=c(5,6,4)

完全要因計画を作成する
引数levels: 変数の水準のベクトル
varNames: 変数の名前> f.design<-gen.factorial(levels=levels.design,varNames=c("A","B","C"))D 最適計画を作成する
引数
data: 候補リスト
nTrials: 最終的な計画の試行回数
evaluateI: 他の基準に加えて，I を評価して，記録するかどうか
criterion: どの基準か (A，D，I)> fract.designD<-optFederov(data=f.design,nTrials=sum(levels.design),evaluateI=TRUE,criterion="D")値
D: 一般化された分散の k 乗根
A: 平均係数の分散
I: X 上の平均予測の分散
Ge: X 上のミニマックス正規化分散
Dea: 近似理論計画のための D 効率の下方限界
design: 計画
rows: 計画の行番号の数値ベクトル> fract.designD$D$[1] 5.012749$A$[1] 0.3586167

I 最適基準I 最適基準は予測値の分散をデータ領域内において平均 $I$[1] 2.71328した値を最小化する最適化基準である．
′f (x) が因子の組み合わせを表す行ベクトルであるとき，$Ge$次の基準を最小化する．
[1] 0.886∫
I=f ′ (x)(X T X)−1 f (x)dx = tr[(X T X)−1 M ](10) $Dea$[1] 0.879

R

この式で，

∫M=

f (x)f ′ (x)dx(11)R

は計画に依存しない情報行列である．(SAS[4] 参照)

5R での実行例

5.1

手段R のパッケージで，一部実施要因計画を生成する AlgDesign を使用する．具体的にはまず，AlgDesign の中にある関数 gen.factorial を使用して完全要因計画 () を生成する．その後，関数 optFederov() に結果を渡し，A，D，I最適計画を作成する．(Web[5] 参照)$design$
A B C
1-2 -5 -3
52 -5 -3
6-2 -3 -3
10
2 -3 -325
2 3 -326 -2 5 -330
2 5 -356 -2 5 -1
91 -2 -5 395
2 -5 396 -2 -3 3
100 2 -3 3
116 -2 5 3
119 1 5 3120 2 5 3$rows$
[1]
191

595

6 10 25 26 30
96 100 116 119 120

56

6シミュレーション

6.1

解析方法

メータを指定し，エラーの数が繰り返しの回数の 10％に満たない場合のみ結果として採用する．

３つの水準でのフルの行列を作り，その行列をランダムに並び替え，指定した計画の数と基準で最適計画を作成する．作成した最適計画の応答曲面を求め，その推定値を取り出し，指定した繰り返し回数での平均と標準偏差を求める．そして，３つ (１次項、２次項、積項) の標準偏差の平均を総合の値として比較する．総合の値を比較するにあたってそれぞれの項に掛ける係数は，１次項が２次項，積項の影響を受けないように 1，２次項，積項は順に 0.5，0.1 とする．水準数は，応答曲面で２次項まで考慮することと実際的な利用から，３水準から５水準の範囲で検討した．なお，応答曲面を計算する際，２水準しか登場しない場合など全てを推定できない場合があることから，計算を行った後，次元数をチェックし，正しい次元に満たない場合はエラーとして，別に記録するようにした．そのエラーの数が繰り返しの回数の 10 ％に満たない場合のみ標準偏差を採用した．そして，２つ以上の基準でエラーが起こらない範囲で，水準と計画の数を変化させ，繰り返し回数を 100 回として解析を行った．

6.2解析手順

プログラムの説明を手順で示す．
手順１
AlgDesign パッケージを呼び出し，使用する変数の初期化をする．
手順２
指定した水準での完全要因計画を作成し，計画の数と同じ数の一様乱数を生成する．その一様乱数を元に計画行列の行をランダムに並び替え，指定した基準の最適計画を作成する．作成した最適計画とその行の数を取り出す．
手順３
最適計画の行の数まで１次項，２次項，積項の計算を繰り返す．さらに，最適計画の行の数までそれぞれの項に係数を掛けたものを足し，さらに平均 0，標準偏差 0.05 の正規乱数のベクトルを誤差項として足す．

y = 1 × x1 + 1 × x2 + 1 × x3+0.5 × x21 + 0.5 × x22 + 0.5 × x23+0.1 × x1 x2 + 0.1 × x1 x3 + 0.1 × x2 x3+rnorm(1, 0, 0.05)

(12)

6.3

３つの水準が全て同じ場合

水準が (3,3,3) の場合，計画の数が 12 以下だと２つ以上の基準でエラーが出る．また，26 以上だと D，A，I の結果が同じになる．計画の数は 3 置きで解析を行った．なお，計画の数の間の取り方として，半分より小さい場合は採用しなかった．(数の間が 5 の場合，3 は採用するが，2 は採用しないで，7 で解析を行った．)
・水準が (3,3,3) の場合
表 1 水準が (3,3,3) での３つの基準の比較 (縦列:計画の数)
DAI
ベスト
13 0.0285 0.0276
0.0268I
16 0.0230 0.0230
0.0225I
19 0.0195 0.0192
0.0196
A
22 0.0168 0.0170
0.0170D
25 0.0167 0.01651 0.01650
I

D が 1 個，A が 1 個，I が 3 個で，水準が (3,3,3) ではI が最もよくなる傾向がある．
３つの水準が全て同じ場合，ベストの数は次のようになる．
表 2 3 つの水準が全て同じ時のベストの数の比較
D A I(3,3,3) 1 1 3(4,4,4) 3 4 3(5,5,5) 5 5 26.4

２つの水準が同じ場合

２つの水準が同じ場合，ベストの数は次のようになる．

表 3 2 つの水準が同じ時のベストの数の比較
D A I(3,3,4) 2 0 3(3,3,5) 2 2 3(4,4,3) 5 0 2(4,4,5) 3 1 4(5,5,3) 4 0 2(5,5,4) 4 1 2

手順４応答曲面を計算し，推定値の次元数を取り出す．もし次元数が正しい場合はそれぞれの項の推定値を取 6.5 ３つの水準が全て異なる場合り出し，それ以外ならばエラーの数として記録する．
水準が (3,4,5) の場合，計画の数が 20 以下だと２つ以手順２〜手順４の作業を指定した回数分繰り返す．
上の基準でエラーが出る．また，60 だと D，A，I の結果手順５が同じになる．計画の数は 5 置きで解析を行った．NA はそれぞれの項の推定値の平均と標準偏差を求める．さ エラーが出て，３つの項の標準偏差の平均が出なかったらに，項ごとの標準偏差の平均を求める．繰り返し ものである．
回数，水準，初期値，基準，計画の数，係数のパラ・水準が (3,4,5) の場合

表 4 水準が (3,4,5) での３つの基準の比較 (縦列:計画の数)
DAI
ベスト
21 0.01271
NA0.01219
I
25 0.00853
NA0.00884
D
30 0.00798
NA0.00793
I
35 0.00747
NA0.00715
I
40 0.00637 0.00834 0.00629
I
45 0.00602 0.00694 0.00620
D50 0.00551 0.00594 0.00552
D55 0.00541 0.00528 0.00550
A59 0.00529 0.00526 0.00527
A

D が 3 個，A が 2 個，I が 4 個で，水準が (3,4,5) ではI が最もよくなる傾向がある．
6.6

考察

水準の与え方と基準の種類との関係は以下のようになる．重み付き平均は，D だと 9/27+20/40+3/9=1.17 が計画 1 回あたりの平均となる．
表 5 水準の与え方とベストの数との比較
DAI
3 つ全て同じ
9108
2 つが同じ20416
3 つ全て異なる324
総計
重み付き平均321.17160.69

281.14

表 5 より，３つの水準が全て同じ場合は A 最適が最もよくなり，２つの水準が同じ場合は D 最適が最もよくなり，３つの水準が全て異なる場合は I 最適が最もよくなることが分かった．また，総計と重み付き平均は D 最適が最も高くなり，全体としては，D 最適が最もよくなるということが分かった．
次に，項ごとの標準偏差の平均の値でのベストの数と基準の種類との関係は以下のようになる．
・３つの水準が全て同じ場合
表 6 3 つの水準が全て同じ時の項ごとのベスト数比較
D A
I
1 次項 11 11 52 次項 11 79
積項999
総計312723

総計は，表 5 では A 最適が最もよくなったが，個別にみた表 6 では D 最適が最もよくなった．これは，A 最適が最もよくはならなくときでも，それに近い値を出していることが要因として考えられる．
項ごとの水準の与え方と基準の種類との関係は 以 下 の よ う に な る ．重 み 付 き 平 均 は ，D だ と31/81+51.5/120+9/27=1.15 となる．

表 7 項ごとの水準の与え方とベストの数との比較
DAI
3 つ全て同じ
3127
23
2 つが同じ51.5 23.5
45
3 つ全て異なる
9612

1 次項
2 次項
積項

31.5
3525

24.5
1121203030

総計
重み付き平均91.5
1.15

56.5
0.75

80
1.10

総計としては，表 5 でみても，表 7 でも D 最適が最もよくなった．よって，応答曲面の観点から考えると D 最適基準が最も適していると考えられる．解析を行う前は，
A 最適が元々応答曲面を意識した計画であることから，A最適が最もよくなるのではないかと考えられたが，実際に解析を行ってみると D 最適が最もよくなり自らの予想とは反した結果になった．
基準 A は，他の２つの基準 D，I と比べて少ない計画の数ではエラーが出やすく，若干劣るということが分かった．解析を行うにあたって，水準を増やすと，比較をするのに最低限必要な計画の数が増えてしまう傾向がある．また計画の数が多いと，誤差は少なくなる傾向がある．よって，水準を増やさないほうが計画の数は少なくて済むが，水準を減らすと誤差が大きくなるという問題がある．また逆に水準が多いと，水準が少ない時と同じくらいの誤差になるのに計画の数が多くなってしまう．(例えば，水準が (3,3,5) の時に計画の数が 19 で済んだのが，水準が(4,4,4) の時に 55 必要になる．) よって，解析を行う人にとっての誤差の許容範囲により，水準と計画の数の決定の仕方が変わってくる．

7

おわりに本研究を通して，計算機支援の実験計画である最適計画についてある程度理解できたと思う．最適計画の方法については，日本語の文献が少なく，英語の文献を読むところから始まったこともあり，理論の学習に多くの時間を費やしてしまったところがあった．また，解析基準の種類にある，G 最適基準と Q 最適基準は計算コストがかかるため取り上げられなかったのは残念である．

参考文献[1] 山田秀: 実験計画法，日科技連，2004．
[2] 轟章: 応答曲面法，http://todoroki.arrow.jp/response/responsesurface.pdf[3] A.C.Atkinson,A.N.Donev: Optimum Experimental
Designs，Oxford Science Publications，1992．
[4] JMP10 実験計画法 (DOE) ガイド，SAS，2012．
[5] AlgDesign R：http://rss.acs.unt.edu/Rdoc/
library/AlgDesign/html/00Index.html