非線形回帰モデルにおける予測情報量規準
九州大学大学院数理学府 金大柱
中央大学理工学部 小西貞則

1

はじめに予測分布とは，将来新たに観測されるデータの予測を行うことを目的として，現在得られているデータとモデルのパラメータに関する事前分布に基づき，ベイズアプローチにより導出される将来新たに観測されるデータに関する分布 (以下，ベイズ型予測分布) である．予測分布は，回帰分析や判別分析など様々なモデルへの応用が試みられている．
Kitagawa (1997) では，分散共分散行列を既知とした下での正規線形モデルに対し，係数に自然共役な事前分布である正規分布を仮定することによって，ベイズ型予測分布が正規分布で与えられることを示した．さらに，情報量規準の観点に基づき，予測分布に基づくモデル評価基準の導出を行い，それを予測情報量規準 (Predictive Information Criterion; PIC) とした．しかし，モデルの分散を既知と仮定しているため理論的には予測分布や PIC の導出が正確に行われているが，一般には，モデルの分散が既知という場面は少なく，実用的な面での問題点が残っていると言える．
本研究では，基底展開法に基づく非線形回帰モデルを仮定し，予測情報量規準の導出を行った．
誤差分散を未知パラメータとし，事前分布として逆ガンマ分布を仮定した場合，予測分布が t 分布となる (Denison et al. (2002), 中妻 (2003)) ため，PIC の導出は困難となるが，ラプラス近似(Tierneay and Kadane (1986), Davison (1986), Kass and Raftery (1995)) を用いることによってPIC を解析的に導出し，基底関数の個数と調整パラメータを選択する手法について考察した．また，数値実験を通して，他のモデル評価基準との比較を行った．

2

基底展開法に基づく非線形回帰モデルにおける予測分布事後分布の導出の際には，自然共役な事前分布を用いることにより，事後分布も事前分布と同じ形で解析的に表現することが出来る (小西・北川 (2004))．このような性質を用い，Bishop (2006)
は，線形回帰モデルの分散を既知とした下で，回帰係数の事前分布として正規分布を仮定し事後分布の導出を解析的に行っている．さらに，分散共分散行列を既知とした線形回帰モデルに対し，このような事前分布を仮定すると，予測分布も解析的に求められることができ，正規分布となることを示した．
線形回帰モデルの回帰係数と分散を，それぞれ，未知のパラメータとして扱う場合の予測分布の導出は，Denison et al. (2002), 中妻 (2003) などで研究されている．このような場合には，回帰係数の事前分布に対して正規分布を仮定し，また，分散の事前分布として逆ガンマ分布を仮定することにより，これらの事前分布は自然共役な事前分布となっており，事後分布も，それぞれ，正規分布と逆ガンマ分布で与えられることが分かる．さらに，事前分布を上記のように仮定した場合のベイズ予測分布の導出を行い，予測分布がスチューデントの t 分布として，解析的に求められることを示した．

12.1

基底展開法に基づく非線形回帰モデル

1 次元目的変数 Y と p 次元説明変数ベクトル X = (X1 , X2 , · · · , Xp )T に関して観測された n 組のデータを {(xi , yi ); i = 1, 2, · · · , n} とする．xi に対して観測された yi は，真の値 u(xi ) に誤差変数 εi が加わって
yi = u(xi ) + εi ,i = 1, 2, · · · , n(1)で与えられたと仮定する．そこで，一般的に回帰モデルとは観測された n 組のデータ {(xi , yi , i = 1, 2, · · · , n)}に対して，

i = 1, 2, · · · , nyi = u(xi ; w) + εi ,(2)で与えられるものとする．この真の構造 u(x) を m 次元パラメータベクトル w = (w1 , w2 , · · · , wm )Tによって特徴づけられるモデル u(x; w) で近似することを考える．
基底展開法に基づく非線形回帰モデルは，回帰関数として，基底関数と呼ばれる既知の非線形関数 bj (x) (j = 1, 2, · · · , m) とパラメータベクトルの線形結合，u(x; w) = w0 + wT b(x)(3)を考える．ここで，b(x) = (b1 (x), b2 (x), · · · , bm (x))T は m 次元基底関数ベクトルとし，w =(w1 , w2 , · · · , wm )T は未知の m 次元パラメータベクトルとする．ここで，w0 は切片である．

y i = w0 +m∑

wj bj (xi ) + εi = w0 + wT b(xi ) + εi ,i = 1, 2, · · · , n(4)j=1

まず，誤差変数 εi が互いに独立に平均 0，分散 σ 2 の正規分布 N (0, σ 2 ) に従うとすると，基底展開法に基づく非線形回帰モデルは，
 {f (yi |w, σ 2 ) =11(2πσ 2 ) 2exp −}2 yi − wT b(xi )2σ 2(5)で表わされる．ここで，推定すべきパラメータは w と σ 2 である．基底展開法に基づく非線形回帰モデルは，基底関数行列 B を用い以下のように表すことができる．

y|w, σ 2 ∼ N (Bw, σ 2 I).

2.2(6)

パラメータの事前分布

本研究では，回帰係数パラメータと分散パラメータの事前分布を，それぞれ，以下のようハイパーパラメータを持つ，正規分布，逆ガンマ分布とする．

w|σ 2 ∼ N (w0 , σ 2 A−1 )()ν0 λ 0σ 2 ∼ IG,.
2 22(7)(8)2.3

事後分布の導出

基底展開法に基づく非線形回帰モデルに対し，回帰係数と分散のそれぞれのパラメータの事前分布が与えられた時，事後分布は次のように定義される．

f (y|w, σ 2 )π(w, σ 2 ).
f (y|w, σ 2 )π(w, σ 2 )dwdσ 2π(w, σ 2 |y) = ∫ ∫(9)

また，尤度関数とパラメータの事前分布が (6)，(8) 式のように与えられた時，パラメータの事後分布は次のような計算により導出することが出来る．
(9) 式の分子は，f (y|w, σ 2 )π(w, σ 2 )= (2π)−m−n22[× exp −(ν2 −m−n−( 20 +1)22(σ )|A−1 − 12|λ02) ν02(Γ ν20)1 {λ0 + (y − B ŵMLE )T (y − B ŵMLE )2σ 2+(w0 − ŵMLE )T{× exp −{−1A−1T+ (B B)

1(w − ŵn )T Â−1n (w − ŵn )2σ 2}−1}](w0 − ŵMLE )}(10)として計算することが出来る．また，事後分布の分母は次のように表すことができる．
∫ ∫f (y|w, σ 2 )π(w, σ 2 )dwdσ 2(= (2π)−n2−1 − 12|A|) ν0λ02(2(Γ ν2012Γ) |Ân | × (λ̂n2ν̂n2)) ν̂n .

(11)

2

ここで，{λ̂n = λ0 + (y − B ŵMLE )T (y − B ŵMLE ) + (w0 − ŵMLE )T A−1 + (B T B)−1}−1(w0 − ŵMLE )(12)とした．よって，事後分布は次のように計算することが出来る．

f (y|w, σ 2 )π(w, σ 2 )f (y|w, σ 2 )π(w, σ 2 )dwdσ 2[]
1−m2 −1 − 12
T2−12= (2π) |σ A | exp − (w − ŵn ) (σ Ân ) (w − ŵn )2π(w, σ 2 |y) =∫∫(×λ̂n2(Γ) ν̂n2( n+ν02 −ν̂n2) (σ )2)+1[]λ̂nexp − 2 .
2σ(13)

これより，回帰係数と分散の事後分布は，それぞれ，以下のようなパラメータを持つ，正規分布，逆ガンマ分布で与えられることが分かる．
(ν̂n λ̂nw|σ , y ∼ N (ŵn , σ Ân ), σ |y ∼ IG,2 22223)

.

(14)

予測分布の導出

2.4

予測分布は次のように計算することが出来る．
∫ ∫f (z|w, σ 2 )π(w, σ 2 )dwdσ 2f (z|y) =∫ {}f (z|w, σ 2 )π(w|σ 2 , y) π(σ 2 |y)dσ 2=∫f (z|σ 2 , y)π(σ 2 |y)dσ 2 .

=(15)予測分布 f (z|y) を求めるために，まず，f (z|σ 2 , y) を求める．

−n−m222f (z|σ , y) = (2π)

2 −n−m22(σ )12|Ân | × (2π)

m2{2σ (B[T}−1 21

B + Â−1n )]{}−11× exp − (z − B ŵn )T σ 2 (B Ân B T + I)
(z − B ŵn )2(16)次に予測分布 f (z|y) を求める．予測分布 f (z|y) は (15) 式より，以下のように求めることができる．
∫f (z|y) =f (z|σ 2 , y)π(σ 2 |y)dσ 2){ (=()nν̂nn + ν̂n/Γ
(πν̂n ) 2
Γ221× 1 + (z − B ŵn )Tν̂n{}λ̂n×(B Ân B T + In )ν̂n

−( n+ν̂n )}−1λ̂n(B Ân B T + I)ν̂n− 122(z − B ŵn )

.

(17)

(17) 式より，予測分布 f (z|y) は以下のようなパラメータを持つステューデントの t 分布に従うことがわかる．
z|y ∼ St(µ̂, Σ̂∗ , ν̂n )

(18)

ここで，

µ̂ = B ŵnλ̂n(B Ân B T + In )ν̂n= (B T B + A)−1Σ̂∗ =Ânν̂n = n + ν0{λ̂n = λ0 + (y − B ŵMLE )T (y − B ŵMLE ) + (w0 − ŵMLE )T A−1 + (B T B)−1}−1(w0 − ŵMLE )ŵn = (B T B + A−1 )−1 (B T B ŵMLE + Aw0 )ŵMLE = (B T B)−1 B T y(19)とする．

3

予測分布に対するラプラス近似

ラプラス近似の目的は，連続変数の集合上に定義される確率密度分布に対しガウス分布による近似を見出すことである．ベイズ予測分布は以下のように表すことができる．
∫h(z|y) ={}exp[n q(w, σ 2 |y) + n−1 log f (z|w, σ 2 ) ]dwdσ 2∫.
exp {nq(w, σ 2 |y)} dwdσ 24(20)小西・北川 (2004) によると，予測分布に対するラプラス近似を用いることにより，上式の q(w, σ 2 |y)のモードを用いることにより，予測分布の正規分布による近似を得ることができる．q(w, σ 2 |y) のw と σ 2 に関するモードをそれぞれ

ŵmode = arg max q(w, σ 2 |y),2σ̂mode= arg max q(w, σ 2 |y)w(21)σ2とする．ここで(λ02) ν02ν0nmn1m1log (2π)− 2 − 2 (σ 2 )− 2 − 2 −( 2 +1) |A−1 |− 2 ( ν0 )nΓ 2q(w, σ |y) =2[× exp −1 {λ0 + (y − B ŵMLE )T (y − B ŵMLE )2σ 2{+(w0 − ŵMLE )T A−1 + (B T B)−1{}−1}](w0 − ŵMLE )}]1(w − ŵn )T Â−1n (w − ŵn )2σ 2(22)arg max q(w, σ 2 |y) = ŵn = (B T B + A−1 )−1 (B T B ŵMLE + Aw0 )(23)

× exp −となる．したがって，w

次に，2σ̂mode= arg max q(w, σ 2 |y)(24)σ2を求めるため，q(w, σ 2 |y) の σ 2 に関する 1 階微分を計算する．σ 2 に依存しない項を C とおくと，()1 (m + n + ν0 )+ 1 log σ 2 + Cn2[11 {+ − 2 λ0 + (y − B ŵMLE )T (y − B ŵMLE )n2σq(w, σ 2 |y) = −{+(w0 − ŵMLE )T A−1 + (B T B)−1}−1}](w0 − ŵMLE ) + (w − ŵn )T Â−1n (w − ŵn )(25)と計算することができる．よって，(25) 式の σ 2 に関する 1 階微分を計算すると，∂q(w, σ 2 |y)∂σ 2

1= −n{}{1(m + n + ν0 )1+1+λ0 + (y − B ŵMLE )T (y − B ŵMLE )2σ 2 2n(σ 2 )2{+(w0 − ŵMLE )T A−1 + (B T B)−1}−1となる．したがって，q(w, σ 2 |y) の σ 2 に対するモードは以下のように計算することができる．

∂q(w, σ 2 |y)=0∂σ 22⇔ σ̂mode={
1λ0 + (y − B ŵMLE )T (y − B ŵMLE )m + n + ν0 + 2+(w0 − ŵMLE )T{−1AT−1+ (B B)}−1}(w0 − ŵMLE ) + (w − ŵn )TÂ−1n (w − ŵn )

.
(26)

5}

(w0 − ŵMLE ) + (w − ŵn )T Â−1n (w − ŵn )これは，以下の式と同等である．

arg max q(w, σ 2 |y) =σ2{}1
λ0 + (ŵn − w0 )T A(ŵn − w0 ) + (y − B ŵn )T (y − B ŵn ) .
m + n + ν0 + 2

よって，ラプラス近似を用いることにより，予測分布は，
22h(z|y, ŵmode , σ̂mode) = f (z|ŵmode , σ̂mode)(1 + Op (n−1 ))=2N (B ŵmode , σ̂modeI)(1 + Op (n−1 ))(27)(28)と近似することが出来る．

基底展開法に基づく非線形回帰モデルにおける予測情報量規準 PIC の4導出
ベイズアプローチに基づく予測情報量規準 PIC (Kitagawa, 1997) は予測分布とデータを生成した真の分布との K-L 情報量の推定量として導かれた．そのとき，予測情報量規準 PIC は，以下のように定義される (Kitagawa (1997))．
22PIC = −2 log h(z|y, ŵmode , σ̂mode) + Bp (q(·), ŵmode , σ̂mode)

4.1(29)バイアス補正項

2
ベイズ型予測分布に対する対数尤度関数 log h(z|y, ŵmode , σ̂mode) は，

11n22I| − 2log h(z|y, ŵmode , σ̂mode) = − log(2π) − log |σ̂mode(z − B ŵmode )T (z − B ŵmode )(30)222σ̂mode2となり，ベイズ型予測分布 h(z|y, ŵmode , σ̂mode) に対する対数尤度と平均対数尤度の差の期待値は，以下のように計算できる．
2Bp (q(·), ŵmode , σ̂mode)[[]]

22= Eq(y) log h(y|y, ŵmode , σ̂mode) − Eq(z) log h(z|y, ŵmode , σ̂mode)= −122σ̂mode{[[tr Eq(y) (y − B ŵmode )(y − B ŵmode )T − Eq(z) (z − B ŵmode )(z − B ŵmode )T]]}

(31)

4.2

真の分布に対する仮定

バイアス補正項を導出するため，真の分布に対し，次のような仮定を置く．真の分布 q(z) があるw∗ ∈ Rm ，σ 2∗ ∈ R1 によって q(z) = f (z|w∗ , σ 2∗ ) で与えられると仮定する．この時，z|w∗ , σ 2∗ ∼N (Bw∗ , σ 2∗ I) とする．

6

4.3

予測情報量規準 PIC

また，パラメータの事前分布の分散共分散行列は A = (nλ)Im とする．この時，ベイズ予測分布に対する平均対数尤度は以下のように計算出来る．
{

Eq(z) (z − B ŵmode )(z − B ŵmode )T}{= Ef (z|w∗ ,σ2∗ ) (z − Bw∗ )(z − Bw∗ )T}+(Bw∗ − B ŵmode )(Bw∗ − B ŵmode )T(32)

ここで，表記 △w ≡ w∗ − w0 ，α = nλ とすると，Bw∗ − B ŵmode = (BB T + αIn )−1 BB T (Bw∗ − y) + α(BB T + αIn )−1 B△w(33)

y − B ŵmode = α(B T B + αIm )−1 (y − Bw∗ ) + α(B T B + αIm )−1 B△wと計算することができる．したがって，バイアス補正項は以下のように計算することができる．
[[]

Eq(y) Eq(z) (z − B ŵmode )(z − B ŵmode )T − (y − B ŵmode )(y − B ŵmode )T[[]]= Ef (y|w∗ ,σ2∗ ) Ef (z|w∗ ,σ2∗ ) (z − B ŵmode )(z − B ŵmode )T − (y − B ŵmode )(y − B ŵmode )T= 2σ 2∗ (BB T + αIn )−1 (BB T )](34)

したがって，バイアス補正項の計算は以下のようになる．
[[22Eq(y) log h(y|y, ŵmode , σ̂mode) − Eq(z) log h(z|y, ŵmode , σ̂mode)]][[[]]]
1TTtr
E(y−Bŵ
)(y−Bŵ)−E(z−Bŵ
)(z−Bŵ)
mode
mode
mode
modeq(y)q(z)22σ̂)([]σ 2∗tr B(B T B + nλIm )−1 B T(35)= −2σ̂mode= −よって，ベイズ型予測分布のモデル評価基準である，予測情報量規準 PIC は以下のように計算することが出来る．

2PIC = n log(2π) + n log(σ̂mode)+4.41(y − B ŵmode )T (y − B ŵmode )2σ̂mode()[]σ 2∗T−1 T+2trB(BB+nλI)B.
n2σ̂mode(36)

数値実験

本節では，シミュレーションによって，提案する予測情報量規準 PIC と従来の一般化情報量規準(GIC) とベイズ型一般化情報量規準 (GBIC) および Eilers and Marx (1996) による方法 (mAIC)を，基底関数との個数と平滑化パラメータおよび事前分布のハイパーパラメータの選択を通して比較検討する．PIC に含まれる真の分散 σ 2∗ に関しては，最尤推定量を代入し数値実験を行った．
ここでは，選択した基底関数の個数と平滑化パラメータの平均と標準偏差，推定曲線と将来の新たなデータとの予測二乗誤差に基づく比較を行う．データは，yα = w(xα ) + εα に基づいて発生させ，真の曲線 w(x) に対しては次の関数を仮定した．

w(x) = 1 − 48x + 218x2 − 315x3 + 145x4 .

7(37)x の定義域は [0, 1] とし，データは等間隔に取った．誤差項 εα に対して正規分布 εα ∼ N (0, σ 2 )，σ = 0.2，0.5，1 を仮定し，標本数 n は n = 20，50，100 の 3 つの設定で実験を行った．表 1 ∼ 3 は500 回のシミュレーションに対する結果である．mean は，基底関数の個数と平滑化パラメータの平均，SD はそれらの標準偏差を表し，PSE は予測値 ŷ と新たに観測されたデータ z との予測二乗誤差PSE =n1 ∑{ŷα − zα }2n α=1(38)である．基底関数の個数は 4 ∼ 15 の範囲で調べ，平滑化パラメータの値は，λ = 10−1 から λ = 10−8を 100 等分した λ1 , · · · , λ100 を調べた．
表 1: n = 20

model selection

number of basis
functionmean – SD

hyper parameterof prior distribution
mean – SDpredictive squarederrormean – SD(τ = 0.2)
mAIC
GIC
GBIC
PIC

12.064 – 3.545
12.536 – 3.205
8.996 – 1.396
11.864 – 0.615

404.001 – 1029.419(×10−5 )91.001 – 433.306 (×10−5 )1149.001 – 684.277(×10−5 )3167.001 – 92.672 (×10−5 )1.341 – 0.483
1.361 – 0.486
1.470 – 0.534
2.333 – 0.588(τ = 0.5)
mAIC
GIC
GBIC
PIC

12.044 – 3.711
12.3 – 3.471
8.828 – 0.976
12.992 – 1.548

350.801 – 705.816 (×10−5 )104.401 – 375.210 (×10−5 )1603.401 – 498.503 (×10−5 )2858.001 – 200.390 (×10−5 )8.254 – 2.981
8.423 – 2.476
7.900 – 2.437
7.788 – 2.437(τ = 1)
mAIC
GIC
GBIC
PIC

12.208 – 3.532
12.368 – 3.422
8.688 – 1.250
13.81 – 1.815

299.401 – 565.171 (×10−5 )114.201 – 338.478 (×10−5 )3280.001 – 2596.629 (×10−5 )2357.601 – 278.563 (×10−5 )32.680 – 11.805
33.377 – 12.132
28.620 – 9.182
26.754 – 8.709

8表 2: n = 50

model selection

number of basis
functionmean – SD

hyper parameterof prior distribution
mean – SDpredictive squarederror
mean – SD(τ = 0.2)
mAIC
GIC
GBICPIC

7.852 – 2.723
8.4 – 2.904
9.676 – 1.596
12.112 – 0.469

56.801 –302.120(×10−5 )30.001 – 188.603 (×10−5 )328.801 – 266.970(×10−5 )1802.801 – 34.562 (×10−5 )2.359 – 0.537
2.374 – 0.537
2.493 – 0.541
3.902 – 0.723(τ = 0.5)
mAIC
GIC
GBIC
PIC

7.326 – 3.032
7.826 – 3.218
10.034 – 1.392
12.804 – 1.335

91.201 – 263.859 (×10−5 )62.001 –217.648 (×10−5 )659.201 – 126.520 (×10−5 )1530.401 – 84.440 (×10−5 )14.506 – 3.259
14.642 – 3.286
15.248 – 3.118
15.591 – 3.152(τ = 1)
mAIC
GIC
GBIC
PIC

7.852 – 3.334
7.94 – 3.315
9.34 – 1.593
12.316 – 2.340

158.001 – 289.696 (×10−5 )93.601 – 240.238 (×10−5 )682.401 –244.724 (×10−5 )1091.201 – 125.035 (×10−5 )57.688 – 13.013
58.280 – 13.208
57.462 – 11.763
56.861 – 11.529

表 3: n = 100

model selection

number of basis
functionmean – SD

hyper parameterof prior distribution
mean – SDpredictive squarederror
mean – SD(τ = 0.2)
mAIC
GIC
GBIC
PIC

7.638 – 2.039
8.046 – 2.377
10.146 – 1.334
12.826 – 0.522

5.401 – 41.885(×10−5 )3.201 – 40.904 (×10−5 )158.401 – 102.813(×10−5 )1098.801–10.899 (×10−5 )4.390 – 0.679
4.409 – 0.686
4.555 – 0.697
6.265 – 0.8857(τ = 0.5)
mAIC
GIC
GBIC
PIC

6.934–2.691
7.238 – 2.807
10.438 – 1.512
13.322 – 1.238

34.401– 121.020 (×10−5 )23.401 – 98.039 (×10−5 )337.801 – 75.912 (×10−5 )870.401 – 47.837 (×10−5 )27.299 – 4.239
27.363 – 4.215
28.501 – 4.257
28.955 – 4.316(τ = 1)mAIC
GIC
GBICPIC

6.95 – 2.896
7.134 – 2.916
9.582 – 1.694
11.854 – 2.327

57.601 – 136.233 (×10−5 )42.801 – 118.305 (×10−5 )329.001 – 77.918 (×10−5 )582.001 – 56.590 (×10−5 )108.240 – 16.684
108.543 – 16.726
110.175 – 16.503
109.693 – 16.402

9

参考文献
[1] Bishop，C．M．(2006)． Pattern Recognition and Machine Learning. Springer.
[2] Davison, A. C. (1986). Approximate predictive likelihood. Biometrika, 73, 323–32.
[3] Denison, D. G. T., Holmes, C. C., Mallick, B. K. and Smith A. F. M. (2002). BayesianMethods for Nonlinear Classification and Regression. Wiley.
[4] Eilers, P. and Marx, B. (1996). Flexible smoothing with B-splines and penalties (withdiscussion). Statistical Science, 11, 89–121.
[5] Hastie, T. and Tibshirani, R. (1990). Generalized Additive models. Chapman and Hall.
[6] Kitagawa, G. (1997). Information criteria for the predictive evaluation of bayesian models.
Communications in Statistics-Theory and Methods, 26, 2223–2246.
[7] Konishi, S., Ando, T. and Imoto, S. (2004). Bayesian information criteria and smoothing
parameter selection in radial basis function networks. Biometrika, 91, 27–43.
[8] 小西貞則，北川源四朗 (2004)． 情報量規準．朝倉書店．
[9] Kullback, S. and Leibler, R. A. (1951). On information and suﬃciency. The Annals ofMathematical Statistics, 22, 79–86.
[10] 中妻照雄 (2003). ファイナンスのための MCMC 法によるベイズ分析, 三菱経済研究所[11] Schwarz, G. (1978). Estimating the dimension of a model. The Annals of Statistics, 6,461–464.
[12] Tierney, L and and Kadane, J. B. (1986). Accurate approximations for posterior moments
and marginal densities. Journal of the American Statistical Association, 81, 82–86.

10