Model Averaging without Non-negativeconstraints
奥井亮 & Guido Kuersteiner
京都大学 & Georgetown University
平成 22 年 11 月 6 日
IBIS2010

概要• 線形回帰モデルにモデル平均手法を応用。
• 予測の平均二乗誤差を最小化する形でモデル平均をするための重み付けベクトルを選ぶ。
• 新しい点： モデル平均をする重み付けベクトルの要素が負になってもよい。
• 利点１： 変数の順序にそれほどの制約がいらない。
• 利点２： 最適な重み付けベクトルの明示的な式をだすことができる。
1

設定{yi, xi}: i.i.d. 標本。yi: １変数。 xi: 回帰変数のベクトルであり、可算無限個の変数からなる。 xi = (x1i, x2i, . . . ).
無限個の回帰変数を含む線形回帰モデルを考える。なお分散均一性を仮定する。

yi = µi + ei,µi =∞∑θj xji,j=1E(ei|xi) = 0,2.
E(e2|x)=σii2

例
この設定はどのような経済学的な状況を表現していると考えることができるか。

• ＧＤＰや株価の予測を多くの経済変数から行う。
• 経済成長率を多くの変数から説明する場合。
• シーヴを使ったノンパラメトリック回帰。

3

近似モデル
さて、真のモデルは、無限個の回帰変数が入っているのでこれを推定することはできない。そこで、ある有限個の変数のみを含むモデルを使って、真のモデルを近似し、y の予測をすることを考える。当然、近似モデルは一つではなく、いくつも考えることができる。
ここでは、0 ≤ k1 < k2 < · · · < km < . . . を整数の数列として、m 番目の近似モ
デルは最初の km 個の回帰変数からなっている場合を考える。
各近似モデルはyi =km∑θj xji + ηmi + ei,

j=1∑∞と書け、ηmi = j=km+1 θj xji が近似誤差である。
4

近似モデルの推定最小二乗法を使って、近似モデルの推定を行う。m 番目の近似モデルによる、µ のxi での推定値 (y の xi での予測値) は、次のように書ける。

µ̂mi =km∑

θ̂j xji,j=1

ここで、θ̂j は、Xkm を、(i, j) の要素が xij である n×km の行列で Y = (y1, . . . , yn)′として、Θ̂m = (θ̂1, . . . , θ̂km )′ は)−1′Θ̂m = Xkm Xkm
Xk′ m Y,(と定義できる。
5

モデル選択とモデル平均
先ほど見たとおり、いくつもの近似モデルを立てることができ、それらの違うモデルからもたらされる予測値は違ってくる。つまり、y の予測をしようとしても、その予測値としていくつもの違った値を考えることができる。

• モデル選択： どれか一つのモデルを選んでその結果を使う。モデル選択は AICや BIC といった情報量基準を立てて行われることが多い。多くの文献があるが、優れた教科書も多いのでそれらを参照されたい。
• モデル平均： 多くのモデルから出てきた結果を組み合わせて、予測を行う。
この発表で取り扱う。教科書としては、Claeskens and Hjort (2008) など。
6

モデル平均µ の M 個のモデルからの推定値を組み合わせる。W = (w1, . . . , wM )′ をモデル∑平均のための重み付けベクトルとする。W ′1M = M
m=1 wm = 1 という条件を課す。(1M は、M × 1 の１からなるベクトルである。)µi のモデル平均推定量はµi(W ) =M∑

wmµ̂mim=1∑である。また、Pm = Xkm (Xk′ m Xkm )−1Xk′ m 、P (W ) = M
m=1 Pm とすると、次のようにも書ける。

µ̂(W ) = P (W )Y =M∑

wmPmY.

m=17重み付けベクトルの選び方
これまでの文献において、いろいろな重み付けベクトルの選び方が提唱されている。

• これまでの方法： ベイズ統計を使った方法。情報量基準の大きさを使う方法。
• ここでは、Hansen (2007) にならい、（条件付）平均二乗誤差を最小化する。
µ = (µ1, . . . , µn)′ かつ、X = (x1, . . . , xn) とすると、Rn(W ) = E((µ̂(W ) − µ)′(µ̂(W ) − µ)|X).

8

平均二乗誤差の形′ (I −P )η
Rn(W ) の式は、Hansen (2007) によって導出されている。am = ηmm mかつ、ηm = (ηm1, . . . , ηmn)′ として、

a 1
 a2

An =  a3
. . .
aMa2 a3a2 a3a3 a3... ...
a M aM

...
...
...
...
...

aM

aM 

aM ,...

aM

k 1
 k1

ΓM =  k1
. . .

k1

k1 k1
k2 k2
k2 k3
... ...
k2 k3



. . . k1
. . . k2 
. . . k3 ... ...
. . . kMとすると、Rn(W ) = W ′(An + σ 2ΓM )Wとなる。
9

最適な重み付けベクトル
最適な重み付けベクトルを Rn(W ) を W ′1M = 1 の条件の下で最小化するものであるとする。
ここで、W の要素を非負であるという制限をおかない。

• あるモデルの重みが負であるとは、そのモデルに含まれている変数を取り除くということと解釈できる。この性質によって、変数の順番付けに関する制約を少なくすることができる。
• W の非負制約がないなら、最小化問題は、明示的な解をもち、計算が簡単になり、また最適な重み付けベクトルの解釈が容易になる。
10

Theorem 1.
W ∗ = (1′M (An + σ 2ΓM )−11M )−1(An + σ 2ΓM )−11M

σ2a −aσ 2 + k1 −k2


2 1
22
σ
σ− 2 a1−a2 + 2 a2−a3
σ + k −kσ + k −k
2 13 222σσ
−+a−aa−a
2+ 2 32+ 3 4σσ.
k−kk−k= 3 24 3...

22σσ−+a−aa−aMσ 2+ kM −1
 σ 2+ kM −2−kM −1M −1 M −2
M −kM −1 


σ2

−+1a−aM−1M
σ 2+ k −kM

M −111W ∗ の式からのいくつかの考察• W ∗ の最初と最後の要素は常に正である。
• W ∗ ∈ [0, 1] × [−1, 1]M −2 × [0, 1].
∗ > 0で• (km+1 − km)(am−1 − am) > (km − km−1)(am − am+1) なら Wmある。もしこの条件が全ての m について満たされるなら、W ∗ ∈ [0, 1]M となり、非負制約をおいた最適解と一致する。この条件は、変数の並べ方が適切である、と解釈できる。
∗ は、m + 1 番目のモデルに入っている変数の説明力のみに依存し、それ• Wmよりも高次のモデルに入っている変数には依存しない。
12

最適な重み付けベクトルの推定先ほどみた最適な重み付けベクトルは、未知のパラメーターに依存しているので、データから計算できない。
ここでは、Hansen (2007) にならい、Mallows 基準を最小化することによって、最適な重み付けベクトルを推定することにする。ここで、ē = (ê1, . . . , êM ), ej =Y ′(I − Pj )Y 、かつ、K = (k1, . . . , kM ) と定義する。最小化する基準は

W ′ē′ēW + 2σ 2K ′W,である。
この最小化問題は、ラグランジェ乗数法で簡単に解くことができる。
13

Theorem 2.
)
1 ′ −1 ( 2′′−1−12′′−12σ K − (1M (ē ē) 1M ) (2 + 2σ 1M (ē ē) K)1M
Ŵ = − (ē ē)
2


k2−k1
2σ ê′ ê −ê′ ê
1 12 2
k3−k2
k2 −k122−σ ê′ ê −ê′ ê + σ ê′ ê −ê′ ê
1 12 22 23 3k−kk−k223243


−σ+σ

ê′2 ê2 −ê′3ê3ê′3ê3−ê′4ê4= .
...


k−kk−kM −1M −1M
M −12
−σ 2 ′
+σêM −2 êM −2−ê′M −1êM −1ê′M −1êM −1 −ê′M êM 
k−k

M
M −1−σ 2 ê′+1′
M −1 êM −1 −êM êM14

Ŵ の式からの考察• Ŵ の最初の要素は常に正である。
• 前に W ∗ ∈ [0, 1] × [−1, 1]M −2 × [0, 1] であることを示したが、Ŵ はその集合に入っている保証はない。あとで見るように Ŵ を使うと実際にはうまくいかないのであるが、このことが、その問題の原因になっていると思われる。特に ê′mêm − ê′m+1êm+1 が小さいなら、対応する Ŵ の要素が非常に大きくなってしまう可能性がある。
• (km+1 − km)(ê′m−1êm−1 − ê′mêm) > (km − km−1)(ê′mêm − ê′m+1êm+1)
なら、m 番目の要素が正になる。もし、この条件が全ての m でみたされるの15なら、非負制約をおいた最適な重み付けベクトルの推定値と、おかないものとは一致する。

• また、E(Ŵ |X) ̸= W ∗ である。ただ、目的関数は不偏である。
• m 番目の要素は、m + 1 番目までのモデルの残差にのみ依存し、それより高次のモデルの残差には依存しない。

y の予測値
y の x での予測値を計算する。d を n × 1 のベクトルで、x = Xk′ d となるものとMする。するとモデル平均による y の予測値は、
ŷ = d′XkM Θ̂ =M∑
m=1(wmd′XkmΘ̂m0)= d′

M∑wmPmY = d′P (W )Y.

m=1となる。

16

y の予測値の表現
P (Ŵ )Yk − k1kM − kM −12P Y
1−σ= σ2 ′ 2PY+1M′′′ê1ê1 − ê2ê2êM −1êM −1 − êM êM+M∑
−1m=2



−σ 2



km+1 − km
km − km−12
 PmY+σ ′′′′êm−1êm−1 − êmêmêmêm − êm+1êm+1P̃m を P̃1 = P1, P̃2 = P2 − P1, . . . , P̃M = PM − PM −1 として定義する。もし、回帰変数が互いに直行しているなら、P̃m は、m 番目のモデルではじめて登場した変数への射影行列である。すると
ê′m−1êm−1 − ê′mêm = Y ′(I − Pm−1)Y − Y ′(I − Pm)Y = Y ′P̃mYとかけるので、M (∑)

km − km−1 ′ŷ = d′P (Ŵ )Y = d′P̃1Y +d P̃mY.
1 − σ2′Y P̃mY
m=217

同等な方法Mallows 基準を最小化する重み付けベクトルを使ったモデル平均推定は、次の方法と同じである。
1. まず、m 番目のモデルに入っていて m − 1 番目のモデルに入っていない変数を、m − 1 番目のモデルに入っている変数に回帰して、その残差をとることによって、変数の直交化を行う。

2. X̃m を、先ほど直交化した、m 番目のモデルに入っていて m−1 番目のモデルに入っていない変数の行列とし、P̃m をその射影行列とする。y をそれぞれの X̃m′ X̃ )−1 X̃ ′ Yに回帰する。すると、各変数の組に対応する係数推定値 β̂m = (X̃mmmを得る。
18

3. β̂m を)(km − km−1∗2β̂mβ̂m = 1 − σ′Y P̃mYとして、縮約する。James-Stein 推定量とは少し違うが非常によく似ている。
(km − km−1 − 2 を使わず、km − km−1 を使っている。)′ )′ での y の予測値を4. X̃ = (X̃1′ , . . . , X̃Mŷ =M∑m=1として計算する。

′ β̂ ∗X̃mm

安定化
重み付けベクトルに非負制約を課さないモデル平均による推定量の性質を、安定的なものにするために、いくつかの方法を考えることができる。

• W ∈ [0, 1] × [−1, 1]M −2 × [0, 1] という制約をおいて、Mallows 基準の最小化を行う。

19

安定化法その２Ŵ を計算をするときに少しの perturbation を入れることによって、Ŵ の要素が非常に大きい値をとることがないようにする。つまり、ある小さな ϵ > 0 を使って、
k2−k1′ê1 ê1 −ê′2ê2+ϵk3 −k222 −k1−σ 2 ê′ ê k−ê′ ê +ϵ + σ ê′ ê −ê′ ê +ϵ1 12 22 23 3σ2...
Ŵϵ = kM −kM −1kM −1 −kM −12−σ 2+σ′′′′êM −2êM −2−êM −1 êM −1+ϵêM −1êM −1 −êM êM +ϵ k−kM
M −12−σ ê′′ ê +ϵ + 1ê−ê
M−1M −1M Mとする。これは、ΥM を (i, j) の要素が M − max(i, j) である M × M 行列として、

W ′(ē′ē + ϵΥM )W + σ 2K ′W,を最小化することで得ることができる。
20

安定化法その３
ここでは、y の予測値を考える。それをŷ + =M∑′ β̂ +X̃mmm=1として計算する。ここで、()km − km−1+2β̂m = max 1 − σ, 0 β̂mY ′P̃mYである。

21

この方法に対応する重み付けベクトルの形では、次のようになる。
()

k2 − k121−σ,0 ,Y ′P̃2Y()()
km+1 − km
km − km−1+22ŵm = max 1 − σ, 0 − max 1 − σ,0 ,′′Y P̃mYY P̃m+1Y()
k −k+ŵM= max 1 − σ 2 M ′ M −1 , 0 .
Y P̃M Yŵ1+ = 1 − max++ ′なお、Ŵ + = (ŵ1, . . . , ŵM) ∈ [0, 1] × [−1, 1]M −2 × [0, 1] である。

モンテカルロ実験
これまで見てきたモデル平均推定量の性質をモンテカルロ実験でみる。
まず、データは次のように生成する。

yi =M
∑

θj xji + ei,j=1

i = 1, . . . , N であり、x1i = 1 で他の回帰変数は xji ∼ i.i.d.N (0, 1) とする。誤差項は ei ∼ i.i.d.N (0, 1) から生成し、回帰変数とは独立である。

22

係数
係数については次の２つの特定化を考える。

• A:

θj = cj −α for j = 1, . . . , M

c は定数である。α がどれぐらい高次の回帰変数が重要であるかを決める。
• B:
θj = 0, for j = 1, . . . , M/2,
θj = cj −α for j = M/2 + 1, . . . , M,
最初に、予測には役にたたない変数が入っていることに注意。つまり、変数の順序付けが間違っている。
23

推定量

• “P” 重み付けベクトルのすべての要素を正であるという制約をおく。
“U” 重み付けベクトルに制約をおかない。
“C” （安定化法１） W ∈ [0, 1] × [−1, 1]M −2 × [0, 1] という制約をおく。
“R” （安定化法 2） W ′(ē′ē+ϵΥ)W +σ 2K ′W を最小化する。なお、ϵ = σ̂ 2とおく。
++“N” （安定化法 3） βmあるいは wmを使う。
• 推定量の性質は、mean of the average squared error (MASE) で評価する:
2
NM
∑

1 ∑ ŷi −θj xji  .
E

N i=1

j=124

結果

• “U” はまったくうまくいかない。
• “R” は “U” よりうまくいく。“R” は ϵ の選び方が重要になってくる。
• “C” はうまくいっている。もし、回帰変数がすべて重要な場合 “P” よりも少し悪くなるものの、M が比較的小さく、R2 が中間的な大きさのときで、回帰変数の並べ方が適切でないときには、“C” は “P” よりもうまくいっている。
• “N” はさらにうまくいっている。回帰変数が適切に並べられている場合でも、“P” と同じぐらいの誤差を返す。
25

結論• 回帰変数が多い場合の線形回帰モデルのモデル平均推定を考えた。
• ここでは重み付けベクトルが非負であるという条件をおかない方法を考えた。
これにより、変数の並べ方にそれほど依存しない方法を考えることができる。
つまりモデル平均において、より多いモデルを少ない計算量で組み合わせることができるのである。
• モンテカルロ実験の結果から、“P” あるいは、“N” を実際には使用するべきであるといえる。どちらを使ったほうがよいかは、回帰変数の並べ方が適切であるかによって、変わってくる。
26