c
オペレーションズ・リサーチ統計的方法における最適化問題
安井 清一
最小二乗法が無制約最小化問題であるように,統計学の中には最適化問題が多く含まれている.解析法の 構成方法には
2
通りがあるように思える.一つは,最小二乗法のように,誤差の二乗和の最小化という問題 を解いて推定量を構成する方法である.もう一つは,算術平均で母平均を推定するといったように,まず統 計量を構成しておいて,その最適性でもって解析法を作る方法である.いずれにしても,最適化問題が含ま れている.統計学というと,数理的な印象を与えるかもしれないので,広くデータ解析を意味する目的とし て,統計学を統計的方法と呼ぶ.本稿ではその中に存在する最適化問題について考える.キーワード:最小二乗法,無制約最小化問題,罰則付き最小二乗法,最良線形不偏推定量,等式制 約付き最適化問題,一般化最小二乗法,最適計画
1.
はじめに統計学,統計的方法の目的は,平たく言うと,母集団 からデータを収集して解析することによって,母集団 の集団としての特徴を探ることである.目的はそうで あっても,データ解析はある種の最適化問題を解くこ とで行われる場面が多く存在するし,解析法の性質を 説明するうえでは最適性が重要になってくる.このよ うに最適化問題は統計的方法の重要な手段である.統 計的方法を説明するうえでは最適化問題を前面に押し 出すことはないが,本稿では,最適化問題を意識して 統計的方法を説明してみたい.
まず初めに馴染み深い単回帰モデルにおける最小二 乗法から入り,重回帰モデル,リッジ回帰における最 適化問題へ行く.その後,統計的方法における推定に ついて,推定量の性質に関して最適化問題を考える.
2.
回帰モデルにおける最適化問題2.1
最小二乗法による単回帰分析統計的方法においても,
OR
においても最も基礎的な 事項は単回帰モデルの回帰係数を最小二乗法によって 求めることだろう.本稿ではまず,ここから始めたい.バ ネ の フック の 法 則 の 実 験 や ,オ ー ム 抵 抗 の 実 験 な ど で ,実 験 や 観 察 に よって ペ ア の デ ー タ
( x
1, y
1) , . . . , ( x
n, y
n)
が取られているとしよう.単回 帰分析では,ペアのデータ( x, y )
間にy
i= β
0+ β
1x
i+ ε
i, i = 1 , . . . , n
やすい せいいち 東京理科大学理工学部
〒
278–8510
千葉県野田市山崎2641
という線形関係がある(線形モデル)と仮定して,パ ラメータ
β
0およびβ
1の値をn
個のデータから求める.統計学においては,データからパラメータを求めるこ とを推定という.また,
n
個の誤差ε
i, i = 1 , . . . , n
に 確率分布(通常は正規分布)を仮定し,傾きβ
1があ るか(β
1= 0
でないかどうか)を調べるための検定,線形関係や誤差の仮定が妥当かどうかなどを調べる回 帰診断と呼ばれる方法までを含めて,統計学では回帰 分析と呼んでいる.詳しくは,佐和
[1]
やDraper and Smith [2]
などを参照してほしい.パラメータ
β
0およびβ
1の推定は,誤差の2
乗和を 最小とするβ
0およびβ
1をそれらの推定値β ˆ
0およびβ ˆ
1とする,おなじみの最小二乗法によって行われる.すなわち,
( ˆ β
0, β ˆ
1) = arg min
β0,β1
n i=1( y
i− β
0− β
1x
i)
2であり,無制約最小化問題である.これは凸関数の最 小化なので,
⎧ ⎪
⎪ ⎪
⎨
⎪ ⎪
⎪ ⎩ dL dβ
0= −2
n i=1( y
i− β
0− β
1x
i) dL
dβ
1= − 2
n i=1( y
i− β
0− β
1x
i) x
iより,連立方程式
⎧ ⎪
⎪ ⎪
⎨
⎪ ⎪
⎪ ⎩
ni=1
( y
i− β ˆ
0− β ˆ
1x
i) = 0
ni=1
( y
i− β ˆ
0− β ˆ
1x
i) x
i= 0
の解
( ˆ β
0, β ˆ
1)
が( β
0, β
1)
の推定値である.この連立方程式を正規方程式という.これを解いて,
β ˆ
0= ¯ y − β ˆ
1x , ¯ β ˆ
1= S
xy/S
xx,
ただし,x ¯ =
ni=1
x
i/n, y ¯ =
ni=1
y
i/n, S
xx=
ni=1
( x
i− x ¯ )
2, S
xy=
ni=1
( x
i− x ¯ )( y
i− y ¯ )
である.2.2
重回帰モデル単回帰モデルにおいて,
x
iを説明変数,y
iを目的変 数という.単回帰モデルは説明変数が一つの場合であ るが,通常はいくつか取ることができて,y
i= β
0+ β
1x
i1+ · · · + β
px
ip+ ε
i, i = 1 , . . . , n
を重回帰モデルという.重回帰モデルの場合,説明変 数が多いので行列表記が有用である.目的変数の値か らなるベクトルをy = ( y
1, · · · , y
n)
とし,切片(の 係数“1”
)と説明変数の値からなる行列X =
⎛
⎜ ⎜
⎝
1 x
11· · · x
1p. .. .
.. .
..
1 x
n1· · · x
np⎞
⎟ ⎟
⎠ ,
偏回帰係数からなるベクトル
β = ( β
0, β
1, . . . , β
p)
, 誤差のベクトルε = ( ε
1, . . . , ε
n)
を用意すると,重回 帰モデルはy = Xβ + ε
のように書ける.偏回帰係数の推定値
β ˆ
は,ε
ε
を最 小にするβ
で求める.すなわち,β ˆ = arg min
β
(y − Xβ)
(y − Xβ) (1)
という無制約最小化問題である.これを幾何学的解釈に よって解く方法もあるが,S (β) = (y−Xβ)
(y−Xβ)
をβ
で微分して解く.一般的に∂ ( β
a ) /∂β = a
,およ び,∂ ( β
Aβ ) /∂β = 2 Aβ
であるので,∂S ( β ) /∂β =
−2X
y + 2X
Xβ
となる.したがって,X
X β ˆ = X
y
を満たす
β ˆ
が(1)
の解である.これは正規方程式の行 列表現である.p + 1
次正方行列X
X
が正則である とき,β ˆ
は一意に定まり,β ˆ = (X
X )
−1X
y
である.ところで,データ数
n
が偏回帰係数の数p + 1
よ りも少ないとき,X
X
は正則でない.また,説明変 数間に相関が強いときや,線形関係が隠れているとき,|X
X|
が小さくなり,逆行列が計算上,不安定にな る.計算上でなくても,誤差を確率変数としたとき,β ˆ
の分散が大きくなり,実際にはあまり有用でなくなる 場合もある.また,n = p + 1
ではすべてのデータをy ˆ = ˆ β
0+ ˆ β
1x
1+ · · · + ˆ β
px
pが通るため,すでに手持 ちのデータに対しての当てはまりはよいが,将来,生 じるデータに対する値については、かなり食い違うと いう現象も生じる.このようなことは,ビッグデータのようなさまざま な種類の変数を大量に扱う場面では,しばしば生じる のではないかと思う.多くの異なる種類の説明変数を 利用することに加え,それらの変数の
2
乗項や積項( x
ix
j, i = j )
,ラグなどを説明変数として利用し,さ らに解析のリアルタイム性を重視するとn p
とな ることも十分考えられる.実用に耐えうるモデルを構 築するためには,変数選択などの一工夫が必要である.回帰モデルの活用目的として,第
1
に予測がある.予測とは,将来生じるデータを推測することである.
回帰分析においては,説明変数は定数なので,ある説 明変数に対する将来生じる目的変数の値の推測である.
上述した場面においては,予測の精度が低下するので,
予測の精度に考慮した説明変数の選択(変数選択,モ デル選択)も必要である.したがって,予測を目的に する場合は,予測の精度を最適化する説明変数の集合 を求めるという問題になる.予測の精度にはいくつか の考え方があるが,赤池情報量規準
(AIC)
が有名であ る.また,統計的機械学習の分野では予測の精度を汎 化能力と呼び,一部のデータのみを用いてモデルを構 築した後,残ったデータに当てはめを行うクロス・バ リデーションという方法で汎化能力を最適化するもの が多いように感じる.赤池情報量規準を代表とする情 報量規準に基づいたモデル選択の詳細については小西 と北川[3]
などがある.以上のように,データから求めた重回帰モデルを実 用において活用するために,変数選択は重要である.変 数選択は
p
個の説明変数からそれらの部分集合を求め る組合せ最適化であるが,実際には2
p通りがモデル の候補ではない.モデルには解釈上の妥当性を考慮し,たとえば,二次項
x
2を含むならば一次項x
も含ませ るという階層性を入れる場合が多く,制約付きの組合 せ最適化となる.変数選択によって目的にあったモデルを構築すると いう方法以外に,罰則項を含んだ最小二乗法による方 法がある.次節ではその代表格であるリッジ回帰を説 明する.
2.3
リッジ回帰と罰則付き最小二乗法ここでは多項式回帰について考える.多項式回帰と は,散布図に
p
次多項式を当てはめることである.す なわち,データ( x
i, y
i) , i = 1 , . . . , n
に対して,y
i= β
0+ β
1x
i+ β
2x
2i+ · · · + β
px
pi+ ε
iを当てはめる.最小二乗法で(偏回帰)係数を求める 場合,前節の話より,
n = p + 1
のとき,すべての点( x
i, y
i)
を通る.明らかに,得られた関数は変動が大き く(上下運動が激しく),予測の精度は高くないこと がわかる.もちろん,AIC
などの規準を用いた変数選 択によって,最適な次数を決めることも可能だが,誤 差の2
乗和に罰則項を追加した罰則付き最小二乗法で,得られる関数の変動を適度に抑えることも可能である.
なお,罰則項は正則化項とも呼ばれ,そのとき,罰則 付き最小二乗法は正則化最小二乗法と呼ばれる.
偏回帰係数の二乗和
pj=1
β
2jを罰則項としよう.定 性的に考えると,より高次の項までモデルに取り入れ るということは0
でないβ
jをより多く含ませるとい うことだから,罰則項は大きくなる.実際にはすべて の係数は非ゼロの値をもつが,大きさが調節されてモ デル全体として適度な変動になると期待される.定式化するためにモデルを少し変形する.切片であ る
β
0は,多項式のy
方向の位置を決めるだけであり,関数の変動とは関係ない.
β
0をなくしたモデルにする ために,各y
iからy ¯
を引き,さらに,各説明変数x
jiの平均値
x ¯
j=
ni=1
x
ji/n
を各説明変数x
jiにおいて引 く.加えて,説明変数x
jiのばらつきは,偏回帰係数の 推定値に影響を与えるので,s
j=
ni=1
( x
ji− x ¯
j)
2 で割っておく.すなわち,y
i− y ¯ = β
1x
1i− x ¯
1s
1+ · · · + β
px
pi− x ¯
ps
p+ ε
iというモデルからスタートする.
z
ij=
xji−¯xjsj として,
モデルの行列表記を新たに
y = ( y
1− y, . . . , y ¯
n− y ¯ )
, β = ( β
1, . . . , β
p)
,
Z =
⎛
⎜ ⎜
⎝
z
11· · · z
1p. ..
z
n1· · · z
np⎞
⎟ ⎟
⎠
と再定義すると,
y = Zβ + ε
である.よって,罰則の大きさを調節するパラメータ
λ > 0
を導入して,S (β) = (y − Zβ)
(y − Zβ) + λβ
β
を最小にするβ
が求めるべき値である.重回帰モデルの ときと同様にして,∂S (β) /∂β = −2Z
y + 2Z
Zβ + 2 λβ
なので,β ˆ = arg min
βS (β)
とすると,(Z
Z + λI) ˆ β = Z
y
の解が求めるべき偏回帰係数の値である.ここで,
I
はp
次の単位行列である.以上より,β ˆ = ( Z
Z + λI )
−1Z
y
となり,これをリッジ回帰推定量という.リッジ回帰は数理計画的な解釈ができる.罰則項は
(y − Zβ)
(y − Zβ)
の最小化に対して,偏回帰係数 の値を束縛する働きがあることから,t > 0
を用いて,min
β( y − Zβ )
( y − Zβ ) s . t . β
β ≤ t
であると考えられる.
λ
を大きくすると罰則が大きく なるため,偏回帰係数β
は最適値として小さい値をと るようになるので,λ
を大きくすることはt
を小さくす ることに対応している.この問題を解くにあたって実 際は,制約条件がアクティブであると仮定される.そ れは,推定においてはt
を小さいほうから少しずつ動 かしてβ
を推定する,クロス・バリデーションなどの 別の基準でアクティブでなくなる前にt
の動きが止ま る,といった理由からであると思われる.結局は等式 制約β
β = t
の下の最小化問題を考えることになる.よって,ラグランジュ関数は
L ( β, θ ) = ( y − Zβ )
( y − Zβ ) + θ ( β
β − t ) , θ > 0
であり,∂L
∂ β = −2Z
y + 2Z
Zβ + 2 θβ = 0
∂L
∂θ = β
β − t = 0
を解けばよい.第
1
行目の方程式から求められるリッ ジ回帰推定量β ˆ = ( Z
Z + θI )
−1Z
y
を第2
式へ代入 し,y
Z(Z
Z + θI)
−2Z
y = t
を得る.よって,t
に 対するリッジ回帰推定量が得られるので,リッジ回帰 における罰則付き最小二乗法は上述の数理計画によっ ても解釈される.2.4
リッジ回帰の一般化 リッジ回帰の罰則項はpj=1
β
j2であった.2
乗和で あると最小化問題の解を解析的に求めることができる が,罰則項としての役割だけを考えると必ずしも2
乗 和である必要がない.そこで,罰則項をpj=1
|β
j|
とし たLeast Absolute Shrinkage and Selection Operator (Lasso)
が提案された.すなわち,Lasso
では,β ˆ = arg min
β
= (y − Zβ)
(y − Zβ) + λ
p j=1|β
j|
を求めることになる.
λ = 0
で通常の最小二乗法にな るので,λ
を大きいほうから0
に向かって少しずつ動 かしながらβ ˆ
を求めるのだが,あるλ
に対して正確に0
になるβ ˆ
j がいくつかあるところがLasso
の特徴で ある.すなわち,リッジ回帰のように罰則付き最小二 乗法を行っているわけだが,同時に変数選択も行って いるという構造になっている.λ
を動かすと,正確に0
になるβ ˆ
jも変化するのだが,どのλ
で止めるかは,クロス・バリデーションで予測の精度を最適化するな どの方法で決められる.
さらに,罰則項を
pj=1
|β
j|
αのように一般化した罰 則付き最小二乗法が考えられている.このときの数理 計画表現は,min
β( y − Zβ )
( y − Zβ ) s . t .
p j=1|β
j|
α≤ t
であり,リッジ回帰の一般化である.この辺りの解説 は
Hastie et al. [4]
にまとめられている.3.
推定に関する最適化問題「ある製品の重さを
n
回測定したら,x
1, . . . , x
n[g]
のようなデータが得られた.その製品の真の重さ
μ [g]
はいくらか?」という問いに対して,真の重さを
n
個 のデータの平均値で求めようとする人がほとんどだろ う.真の重さμ
をデータから言い当てることを推定と いい,推定のための式を推定量,実際に計算した値を 推定値という.推定量および推定値をμ ˆ
と書き,今こ こでは,推定量と推定値とを特に区別しなくても不都 合はないから,μ ˆ = ¯ x =
ni=1
x
i/n
である.しかしな ぜ,平均値を選ぶのであろうか.子どもの頃から言わ れている,データを増やせば平均値は精密になる,す なわち,大数の法則が感覚としてある(どこかで刷り 込まれた?)などが理由であるような気がする.本節では,推定量としての平均値へのこだわりについて考 察する.
3.1
不偏推定量製品の真の重さを
μ
としたとき,n
回測定して得ら れたデータをx
i= μ + ε
i, i = 1 , . . . , n (2)
のように考える.ε
1, . . . , ε
nは測定誤差を表しており,互いに独立で同一な分布に従う確率変数であるとする.
また,期待値と分散はそれぞれ
E [ ε
i] = 0
,V [ ε
i] = σ
2 とする.この仮定の中でE [ ε
i] = 0
が特に重要であり,これは測定に偏りがないということを示しており,標 準試料によって校正がきちんと行われているなどの現 実的な意味が含まれている.測定に偏りがあると真の 値を推定するのが困難になる.また,これらの各仮定 は,校正に加え,測定時においても一定の手順(標準 という)が定められており,それに従って熟達した人 が行った結果であるという意味がある.たとえば,測 定に未熟な人が行った場合,
n
回測定している間に上 達して,だんだんと分散が小さくなるといったような ことは起こらない,ということである.さて,
μ
の推定量である¯ x
の期待値を求めてみよう.x ¯ = μ +
ni=1
ε
i/n
であるので,E [¯ x ] = μ + 1 n
n i=1E [ ε
i] = μ
である.推定量の期待値が推定対象に一致する推定量 のことを不偏推定量という.すなわち,
x ¯
は平均的に 推定対象μ
をとる推定量であり,また,大数の法則を 参照すると,n
を大きくする(データを増やす)と推 定対象μ
に近づくという推定量である.μ
の 不 偏 推 定 量 は 他 に も あ る .た と え ば ,x
2, x
4, . . . , x
2n/2のような添え字が偶数のものだけを用 いた平均値x ¯
1=
n/2k=1
x
2k/n/ 2
もそうである.ゆ えに,x ¯
はμ
の不偏推定量だからよい,とは言い切れ ないのである.そこで次に推定量の分散を比較するこ ととなる.3.2
最小分散(最良)線形不偏推定量 平均値x ¯
の分散は,V [¯ x ] = V
ni=1
ε
in
= 1 n
2 n i=1V [ ε
i] = σ
2n
である.また,同様にして,V [¯ x
1] = σ
2/n/ 2
であ る.よって,x ¯
1よりも¯ x
の分散のほうが小さいことがわかる.どちらの推定量も
ni=1
a
ix
iという形をして おり,線形推定量と呼ばれる.x ¯
はすべてのa
iが1 /n
であり,¯ x
1は奇数添え字のa
iを0
,偶数添え字のa
iを1 /n/ 2
としたものである.また,これらは線形推定 量かつ不偏推定量なので,線形不偏推定量と呼ばれる.ここで,線形不偏推定量の中で,最小の分散をもつ推 定量,すなわち,最小の分散を与える
a
i, i = 1 , . . . , n
を求めてみよう. ni=1
a
ix
i が 不 偏 推 定 量 で な け れ ば な ら な い の で,E
ni=1
a
ix
i= μ
でなければならない.す なわち,ni=1
a
i= 1
の下でV
ni=1
a
ix
i を最 小にするa
i, i = 1 , . . . , n
を求めることとなる.V
ni=1
a
ix
i=
ni=1
a
2iσ
2だから,a1,...,an
min
ni=1
a
2is . t .
ni=1
a
i= 1
という等式制約付き最小化問題を解けばよい.形式的 ではあるがラグランジュ未定係数法を適用する.
a = ( a
1, . . . , a
n)
として,ラグランジュ関数をL ( a, λ ) =
ni=1
a
2i− λ (
ni=1
a
i− 1)
とする.よって,∂L (a, λ )
∂a
i= 2 a
i− λ = 0 , i = 1 , . . . , n
∂L (a, λ )
∂λ =
n i=1a
i− 1 = 0
を解けばよい.
2
ni=1
a
i= nλ
,および,ni=1
a
i= 1
よりλ = 2 /n
である.ゆえに,a
i= 1 /n, i = 1 , . . . , n
であり,x ¯
がμ
の線形不偏推定量の中で最小の分散を もつことがわかった.このような推定量を最小分散線 形不偏推定量,もしくは,最良線形不偏推定量という.μ
の推定量として感覚を通じてx ¯
を選んだが,x ¯
は データのモデル(2)
においてμ
に対する最小二乗法で 得られる推定量(最小二乗推定量)になっている.一 般的な重回帰モデルにおいても,どの誤差も期待値お よび分散が0
およびσ
2(等分散)であり,どの二つの 誤差を見ても無相関であるという条件の下で,最小二 乗推定量は最良線形不偏推定量であることがガウス・マルコフの定理で示される.詳しくは佐和
[1]
などを 参照されたい.最後に,
(2)
において,誤差の仮定を一般的にしたモ デルに対して,μ
の最良線形不偏推定量を求めてみる.すなわち,
V [ ε
i] = σ
i2, i = 1 , . . . , n
,Cov [ ε
i, ε
j] =
0 ( i = j )
とする.行列を用いたほうが便利なので,デー タを行列表記すると,x = ( x
1, . . . , x
n)
とすると,x = μ1 + ε, E [ ε ] = 0, V [ ε ] = Σ (3)
となる.ただし,E [ ε ]
は期待値E [ ε
i]
を列ベクトルに 並べたものであり,V [ ε ]
は対角要素に分散V [ ε
i]
,非 対角要素に共分散Cov [ ε
i, ε
j]
を配置した対称行列で ある.よって,分散の性質より,Σ
は正値定符号行列 である(非負値定符号行列の場合もあるが,ここでは 正定性を仮定している).任意の線形推定量はa
x
で ある.不偏性より,E [ a
x ] = a
E [ x ] = μa
1
であ るからa
1 = 1
でなければならない.推定量の分散V [a
x] = a
V [x]a = a
Σa
より,μ
の最良線形不偏 推定量はmin
aa
Σa s . t . a
1 = 1
の解より得られる.ラグランジュ関数は,L ( a, λ ) = a
Σa − λ ( a
1 − 1)
であり,∂L ( a, λ )
∂a = 2Σa − λ1 = 0
∂L (a, λ )
∂λ = a
1 − 1 = 0
を解けばよい.第
1
番目の式に左からΣ
−1をかけると,2 a−λΣ
−11 = 0
が得られ,さらに左から1
をかけて第2
番目の式に代入すると,2−λ1
Σ
−11 = 0
を得る.よ って,λ = 2 /1
Σ
−11
である.これを,2a−λΣ
−11 = 0
に代入すると,a = Σ
−11 1
Σ
−11
である.よって,
μ
の最良線形不偏推定量はμ ˆ = a
x = ( 1
Σ
−11 )
−11
Σ
−1x
である.なお,これも
(2)
に対する最小二乗推定量と 同様に,(x − μ1)
Σ
−1(x − μ1)
を最小にする
μ
であり,(3)
に対する一般化最小二乗 法による推定量(一般化最小二乗推定量)でもある.3.3
再び単回帰分析(最適計画)2.1
節の単回帰分析では,n
個のデータ( x
i, y
i)
が与 えられ,それに線形式を当てはめることが目的であった.ここでは,その逆を考えたい.よい線形式を得る ためにはどのようなデータをとればよいかを求める.
よい線形式というのを,回帰係数の推定量の分散 が小さいものと定義する.単回帰モデルは
p = 1
の 重回帰モデルなので,回帰係数の推定量はp = 1
としたβ ˆ = ( X
X )
−1X
y
である.推定量の分散 は,V [( X
X )
−1X
y ] = σ
2( X
X )
−1 である.推 定量の分散は行列であり,大小比較ができないので,V [ˆ β]
の行列式にして比較することをよく行う.すな わち,|σ
2( X
X )
−1|
を最小にするX
を求める問題で ある.X
の中身は,データ( x
i, y
i)
のx
iからなるの で,|σ
2(X
X)
−1|
を最小にするために,どのx
でy
を得ればよいかという問題となる.この問題はx
を自 由にどの点にでも取ることができなければ成立しない ので,x
は制御できる変数でなければならない.つま り,これはx
のある点で実験を行いy
を観測するとい う場面に相当する.このことから,望ましいx
,すな わち,望ましい実験点の集合X
を決める問題を最適計 画を求める問題という.通常,実験は実験ができる範 囲があったり,また,理論的に考える場合でも,実験 可能領域を定めないと際限がないので,閉区間[−1 , 1]
を実験点の領域とすることが標準的である.なお,こ こでの計画は
“design”
であり,数理計画でいう“pro-
gramming”
とは漢字は同じだが意味は異なる.ここで,
min
X|σ
2( X
X )
−1|
はmax
X|X
X|
と同じであ ることに注意すると,最適計画を求める問題はmax
X|X
X | s . t . ∀i, x
i∈ [−1 , 1]
である.この問題から得られる最適計画は,
X
X
の行列式
Determinant
を目的関数(最適性の基準)とするところから
D -
最適計画 と呼ばれる.目的関数をtr(X
X)
−1,つまり,推定量の分散の和V [ ˆ β
0]+ V [ ˆ β
1]
を最小にする計画を考えるものもあり,この問題で得 られる計画をA-
最適計画 という.この他にも,E -
最適性やI -
最適性など,さまざまな基準がある.さて,
1
次元ではあるが単回帰モデルについてD-
最 適計画を求めよう.以下,話を簡単にするためn
は偶 数という仮定を追加する.|X
X|
はn
ni=1
x
i ni=1
x
i ni=1
x
2i= n
n i=1x
2i−
ni=1
x
i 2= n
n i=1( x
i− x ¯ )
2= nx
I − J n
x
と変形できる.
J
はすべての要素が1
であるn
次正方 行列である.( I − J/n )
はべき等・対称行列であるので,|X
X| = || (I − J/n ) x||
2である.また,(I − J/n )
は1
が張る空間に対する直交補空間上への射影行列な ので,|X
X|
を最大にするためには,x
を( I − J/n )
の列空間上に取ればよい.このことは,最適なx
をx
∗ とすると,1
x
∗= 0
,すなわち,最適なx
にお けるx
j, j = 1 , . . . , n
の平均値は0
であるという ことである.よって,x
∗= ( x
∗1, . . . , x
∗n)
とすると,|X
X| = n
ni=1
( x
∗i)
2かつ|x
∗i| ≤ 1 , ∀i
であるので すべてのi
についてx
∗i= ± 1
となり,1
x
∗= 0
を満 たすものは,n
を偶数としたので,「x
iの半分を−1
に,もう半分を
+1
」にした計画が単回帰モデルにおけるD -
最適計画である.すなわち,実験領域の端で実験を 行うのがよいということである.このことは,1
次項 のみの線形モデルに対して,p ≥ 2
においても基本的 に成り立つ.実験をもっと意識すると,データ数
n
は実験回数,説明変数の数
p
は実験に取り上げる因子の数である.n
とp
に対して最適計画を求める方法もあるが,「X
の第1
列目を1
とし,残りの要素が±1
であり,各列 が直交する行列」H
sを使う方法もある.このような 行列H
sに従う実験は,すべての実験点が超立方体の 頂点であり,D -
最適計画である.すべての実験点が超 立方体の頂点であるということは,どの因子も− 1
お よび1
の2
条件(これを2
水準という)の実験になる ので実験が行いやすい.そこで,H
sタイプのn
次正 方行列,すなわち,「第1
列目が1
,残りの要素が± 1
のみである直交行列」H
を見つけたいのだが,この問 題は簡単ではない.行列H
はアダマール行列と呼ば れ,n = 1 , 2
,それ以上はn
が4
の倍数のときに存在 する.さらに,「行または列の入れ換え」および「行ま たは列の符号の反転」によって同一になるものを同型 とするとき,n = 1 , 2 , 4 , 8 , 12
においては一つしかな く,n = 16
は五つ存在する.たとえば,n = 1 , 2 , 4
に 対しては,(1)
1 1 1 −1
⎛
⎜ ⎜
⎜ ⎜
⎝
1 1 1 1
1 1 − 1 − 1
1 −1 1 −1
1 −1 −1 1
⎞
⎟ ⎟
⎟ ⎟
⎠
n = 1 n = 2 n = 4
である.また,アダマール行列の構成法はいくつかあ り,その一つとして
H
2k=
1 1
1 −1
⊗ H
2k−1, k = 1 , . . .
(
Sylvester
タイプ)がある.
⊗
はクロネッカー積である.この方法ではn = 12
やn = 16
の残りの四つは作れない.しかし,アダマール行列を実験の計画および統計解析に用いた とき,
Sylvester
タイプとそれ以外のタイプとでは,統 計的な性質が異なる.最適計画の問題を含め,実験を効率的に行い,情報 を効果的に得るためには,どのような実験点を取るか
(実験の計画),また,どのように解析を行えばよいか を示した方法を実験計画法という.実験計画法の良書 は多くあるが,山田
[5]
は,基礎的な方法から最適計 画のようなアドバンスな方法まで,内容が豊富である.実験計画法において忘れてはならないのが田口玄一で ある.アダマール行列は直交配列表実験として,田口 以前から使用されていたが,実験用にアダマール行列 を使いやすく書き直す工夫をすることにより,日本に おいては直交配列表実験が専門家以外の人々に普及し,
製造業では一般的な実験方法となった.直交配列表実験 とは実験回数が制約された中から,最大限に情報を得 るための実験計画であり,製造条件の最適化などに役 立ってきた.また,田口は製品の使用環境で生じるであ ろうノイズを製品開発時の実験で発生させ,ノイズの 影響を緩和して,ノイズに強い製品を設計するため手 法であるロバストパラメータ設計
(Robust Parameter
Design)
を開発した.さらに田口は製品開発・技術開発に関わる新しい実験・データ解析技術を多く開発し,
実験計画法や統計的方法にも影響を与えた.それらの 手法群はタグチメソッドと呼ばれており,世界でも知 られている.タグチメソッドの入門書として立林
[6]
, 実験計画法・統計的方法の側面においては,日本では 宮川[7]
,世界ではWu and Hamada [8]
などがある.4.
おわりに統計的方法の背後には最適化理論がある.サポート ベクターマシンのように,統計的機械学習の手法にお いては,凸計画問題が陽の形で現れている.しかし,最 適化理論は解析法の基礎をなすものの,統計的方法は それ以外の部分も,等しく重要であることにも注意し たい.たとえば,得られたデータが観察によるものな のか,実験によるものなのかによって,同じ回帰分析 でも,最終的に導くことができる結果に違いが出てく る.ランダムな順序で行った実験からは取り上げた因 子と特性(データ)との因果関係を統計的に検証でき るが,そうでない観察から得られたデータからは,基 本的に予測しか保証されていない.また,母集団の規 定,サンプリング方法なども解析目的へどれだけ迫れ るかに影響を与える.解析法,その基礎技術である最 適化理論は,統計的方法の中心に位置するが,その他 の要素が適切であってこそ,解析の目的を達成できる ことに,最後に触れておきたい.
参考文献