日本統計学会誌, 第45巻, 第1号, 143頁-170頁

日本統計学会誌
第 45 巻, 第 1 号, 2015 年 9 月
143 頁 ∼ 170 頁
特集
スタイン・パラドクスとベイズ理論
丸山祐造∗
A Review: The Stein Paradox and Bayes Theory
Yuzo Maruyama∗
多変量正規分布やそれを拡張した球面対称分布のパラメータ推定や予測問題において生じるス
タイン・パラドクスに関して概説する．特に，分散既知のもとでの多変量正規分布の平均ベクト
ルの推定，平均ベクトルの推定と予測分布の関係，線形回帰モデルにおけるパラメータ推定を扱
う．推定量や予測分布の決定理論的な良さを考える際に，ベイズ推定量やベイズ予測分布が重要
な役割を果たすことを概観する．
The Stein paradox on estimation and prediciton under multivariate normal distribution and
spherically symmetric distribution is reviewed. In particular, estimation of a multivariate normal mean with known variance, the relationship of the esimation of mean vector and predictive
density and parameter estimation of linear regression model are treated. Bayesian procedures
play an important role as a remedy for the Stein paradox.
キーワード: スタイン・パラドクス，統計的決定理論，許容性
イントロダクション
1.
本稿では，多変量正規分布やそれを拡張した球面対称分布のパラメータ推定や予測問題
において生じるスタイン・パラドクスに関して，いくつかのトピックを紹介する．ただし，
久保川 (2004) による優れた概説が既に存在するため，研究分野としてバランスの取れた概
説を目指すのではなく，筆者の視点で興味深いと思うトピックに絞って紹介していく．特
に推定量や予測分布の決定理論的な良さを考える際に，ベイズ推定量やベイズ予測分布が
重要な役割を果たすことを概観する．
第 2 節では，分散既知のもとでの多変量正規分布の平均ベクトルの推定について概説す
る．2.1 節で，一般化ベイズ推定量が許容的であるための十分条件について考える．事前
分布の密度関数に，原点に関する球面対称性，微分可能性，正則変動性を仮定したもとで，
∗
東京大学空間情報科学研究センター：〒 277-8568 千葉県柏市柏の葉 5-1-5 (E-mail: [email protected]).
144
日本統計学会誌第45巻第1号 2015
Brown (1971) による十分条件と本質的に同等の十分条件を与える．2.2 節では，標本平均
ベクトルを改良するための十分条件について，優調和性を軸にして議論する．第 3 節では，
予測分布と第 2 節のパラメータ推定の関係を概説する．平均ベクトルの推定において標本
平均ベクトルを改良するベイズ推定量を導く事前分布が，予測分布においてもその良さを
引き継ぐことを確認する．第 4 節では，線形回帰モデルにおけるパラメータ推定において
生じるスタイン・パラドクスを扱う．ある特別な事前分布のクラスを考えるとき，球面対
称分布のもとでの一般化ベイズ推定量が，正規分布のもとで計算された一般化ベイズ推定
量に一致するというロバスト性が生じることを示す．さらに一般化ベイズ推定量のクラス
の中に，Stein (1956) の推定量に関連する簡潔な表現を持つ推定量があることを示す．
2.
分散既知のもとでの平均ベクトルの推定
yij (i = 1, . . . , n, j = 1, . . . , d) が互いに独立に正規分布 N (θj , 1) に従うとき，θ1 , . . . , θd
を同時に推定する問題を考える．正規分布のもとでは，標本平均ベクトル
(∑ n
i=1
yi1 /n, . . . ,
∑n
i=1
)T
yid /n
= (x1 , . . . , xd )T = x
(2.1)
が θ = (θ1 , . . . , θd )T ∈ Rd に関する十分統計量であるので，始めから d 変量正規分布
Nd (θ, n−1 Id ) に従う x だけが得られると考えてよい．また本節の議論では分散行列が既知
であることが本質的なので，簡単のため単位行列とする．つまり x ∼ Nd (θ, Id ) が得られ
たときに，平均ベクトル θ を推定する問題を扱う．
推定量 δ(x) の良さを平均二乗誤差 R(θ, δ) = E[kδ(x) − θk2 ] で評価しよう．もし推定
量 δ に対して「任意の θ に対して R(θ, δ∗ ) ≤ R(θ, δ) であり，ある θ に対して強い不等号
が成り立つような δ∗ が存在する」ならば，δ でなく δ∗ を使うべきである．このとき δ を
非許容的であるといい，δ∗ は δ を改良する推定量という．一方，推定量が非許容的でない
（つまり，他のどんな推定量によっても改良されない）ならば，許容的であるという．もち
ろん許容性は良さの基準としては非常に弱いが，それ故推定量に最低限持って欲しい性質
である．
標本平均に対応する推定量 x は，最尤性，最小分散不偏性を持つ自然な推定量なので，
許容的であるのが当然と思われる．実際 d = 1, 2 のときには許容的である．しかし d ≥ 3
のとき，x を改良する推定量が存在する．直感に反するこの現象は，発見者 Stein (1956)
の名を取って「スタイン・パラドクス」と呼ばれる．本稿では，平均ベクトルの次元 d は
3 以上であることを仮定する．
x を改良するための十分条件を求めるにあたり，スタイン等式とよばれる期待値に関す
る恒等式が大きな役割を果たす．z を正規分布 N (µ, σ 2 ) に従う確率変数，g(z) を絶対連続
スタイン・パラドクスとベイズ理論
145
で E[|g 0 (z)|] < ∞ を満たす関数とする．このとき z の密度関数
φ(z; µ, σ 2 ) = (2πσ 2 )−1/2 exp(−(z − µ)2 /{2σ 2 })
が (z − µ)φ(z; µ, σ 2 ) = −σ 2 {d/dz}φ(z; µ, σ 2 ) を満たすので，部分積分より
∫ ∞
∫ ∞
g(z)(z − µ)φ(z; µ, σ 2 )dz = σ 2
g 0 (z)φ(z; µ, σ 2 )dz
E[(z − µ)g(z)] =
−∞
−∞
(2.2)
0
2
= σ E[g (z)]
が成り立つ．この恒等式 E[(z − µ)g(z)] = σ 2 E[g 0 (z)] がスタイン等式として知られてい
る．スタイン等式を用いて，推定量 x + k(x) が x を改良するための十分条件を与えよう．
x + k(x) の平均二乗誤差は
R(θ, x + k(x)) = E[kx + k(x) − θk2 ]
= E[kx − θk2 + kk(x)k2 + 2(x − θ)T k(x)]
(2.3)
= d + E[kk(x)k2 + 2 ∇T k(x)] = d + E[ξ[k(x)]]
と書ける．ただし，∇ = (∂/∂x1 , . . . , ∂/∂xd )T であり，
ξ[k(x)] = kk(x)k2 + 2 ∇T k(x)
(2.4)
である．(2.3) の３番目の等号がスタイン等式による．従って，
ξ[k(x)] ≤ 0 ∀x
(2.5)
ならば，x + k(x) が x を改良する．例えば kJS (x) = −{(d − 2)/kxk2 }x に対して，ξ[kJS ] =
−(d − 2)2 /kxk2 ≤ 0 である．対応する推定量
δJS = x + kJS (x) = (1 − (d − 2)/kxk2 )x
は James and Stein (1961) による James-Stein 推定量として知られている．
ところで δJS は kxk2 > d − 2 のとき，縮小係数
0 < 1 − (d − 2)/kxk2 < 1
により x を 0 方向に縮小している．従って kθk が小さいときに効果的な推定量であると解
釈できる．しかし，kxk2 < d − 2 のときには，係数は負であり，δJS と x の符号が逆にな
る．この欠点を修正した James-Stein positive-part 推定量
δJS+ = max(0, 1 − (d − 2)/kxk2 )x
146
日本統計学会誌第45巻第1号 2015
+
が δJS を改良するため1) ，δJS は非許容的である．ところが δJS
も非許容的なのである．2.1
節の定理 2.1 にあるように，狭義の事前分布に対するベイズ推定量の集合と広義の事前分布
に対する一般化ベイズ推定量の集合の和集合が，全ての許容的な推定量を包含するからで
+
ある．関数として滑らかでない δJS
はこの和集合には入らない．なお狭義の事前分布に対
するベイズ推定量は許容的であるが，広義の事前分布に対する一般化ベイズ推定量は許容
的とは限らない．例えば，x は広義の一様分布に対する一般化ベイズ推定量であるが，既
に述べた通りその許容性は次元 d に依存する．2.1 節では広義の事前分布に対する一般化
ベイズ推定量が許容的であるための十分条件を扱う．
2.1
許容性
全ての許容的な推定量はベイズ推定量または一般化ベイズ推定量として特徴づけられる
ことが，Brown (1971) によって示された．
定理 2.1 (Brown (1971)) 推定量 δ が許容的であれば，ある非負の測度 ν(dθ) が存在
して
∫
mν (x) =
(2π)−d/2 exp(−kx − θk2 /2)ν(dθ)
(2.6)
とおくとき，全ての x について mν (x) < ∞ で，しかもルベーグ測度に関してほとんどい
たるところで
δ = x + ∇ log mν (x)
(2.7)
と表される．
また Brown (1971) は，
（一般化）ベイズ推定量が許容的であるための十分条件を，周辺
尤度 (2.6) に対する条件として与えた．
定理 2.2 (Brown (1971)) 非負の測度 ν(dθ) に対する周辺尤度 mν (x) が以下の二つ
を満たすとする．
A1 k∇ log mν (x)k が有界である．
A2 J を区分的に微分可能で，kxk ≤ 1 において j(x) ≥ 1 を満たす関数 Rd → R の集合
とする．このとき，kxk ≥ i において ji (x) = 0 を満たし，また
∫
k∇ji (x)k2 mν (x)dx = 0
lim
i→∞
(2.8)
Rd
であるような関数列 ji ∈ J が存在する．
1)
証明は Lehmann and Casella (1998) Theorem 5.5.4 や Anderson (2003) Lemma 3.5.2 を参照のこと．
スタイン・パラドクスとベイズ理論
147
このとき ν(dθ) に対する（一般化）ベイズ推定量は許容的である．
さらに Brown (1971) は，測度 ν が原点に関して球面対称であれば，その（球面対称な）
周辺尤度 m̄ν (kxk) = mν (x) に対する条件
∫ ∞
dr
?
=∞
A2
d−1 m̄ (r)
r
ν
1
で，条件 A2 が置き換えられることも示した．
ただし，定理 2.2 はその一般性ゆえ，証明を完全に理解するのは難しい．本稿では，事
前分布に対していくつかの仮定をおき，平易な証明によって十分条件を与える．その仮定
とは，事前分布 ν(dθ) の原点に関する球面対称性，密度関数
ν(dθ) = π(kθk)dθ
(2.9)
の存在，さらに π(r) の正則変動性および微分可能性である．
正則変動性は漸近的に多項式オーダーで振る舞う関数を扱うのに優れた概念であり，
π(tr)
= tα
r→∞ π(r)
lim
∀t>0
(2.10)
を満たす π はインデックス α の正則変動関数とよばれる．rα のほか対数関数を含む
rα {log r}β ，rα {log r}β {log log r}−γ などを，インデックス α の正則変動関数として統一的
に扱うことが出来る．
また π(r) の微分可能性のもとで，(2.6)，(2.7)，(2.9) に基づいて与えられる（一般化）ベ
イズ推定量
∇
δπ = x + ∫
∫
(2π)−d/2 exp(−kx − θk2 /2)π(kθk)dθ
(2π)−d/2 exp(−kx − θk2 /2)π(kθk)dθ
Rd
Rd
(2.11)
に対して，以下のような別表現が与えられる．右辺の分子において θ ∼ Nd (x, Id ) と見な
してスタイン等式 (2.2) を用いると
∫
∇
Rd
∫
exp(−kx − θk2 /2)
π(kθk)dθ
(2π)d/2
Rd
∫
exp(−kx − θk2 /2)
=
{∇π(kθk)}
dθ
(2π)d/2
Rd
exp(−kx − θk2 /2)
π(kθk)dθ =
(2π)d/2
(θ − x)
を得る．(2.12) によって（一般化）ベイズ推定量の別表現
∫
2
d {∇π(kθk)} exp(−kx − θk /2)dθ
= x + Eθ|x [∇ log π(kθk)]
x + R∫
exp(−kx − θk2 /2)π(kθk)dθ
Rd
(2.12)
(2.13)
が得られる．ここで右辺の期待値は，事後分布に関する ∇ log π(kθk) の期待値である．もと
もと (2.12) 及び (2.13) は，Brown and Hwang (1982) で与えられた表現である．Brown and
Hwang (1982) もまた，事前密度関数に対する微分可能性の仮定によって，Brown (1971)
148
日本統計学会誌第45巻第1号 2015
より平易に十分条件を与えるという意図があった．ただし，証明において不等式の評価が
粗いために，π(kθk) = kθk−(d−2) をはじめとして，Brown (1971) の許容性の十分条件を
満たす事前分布が Brown and Hwang (1982) の許容性の十分条件を満たさない．次の定理
では，Brown and Hwang (1982) の方法を踏襲しつつ，また Brown and Hwang (1982) と
Maruyama (2009) の証明に改良を加えて，一般化ベイズ推定量が許容的であるための十分
条件を与える．その証明に（一般化）ベイズ推定量の別表現 (2.13) が重要な鍵となる．
定理 2.3 (Maruyama (2009)) 事前分布の密度関数 π(kθk) において，π(r) が 2 階微
∫1
分可能であり，原点付近で rπ 0 (r)/π(r) が有界， 0 rd−2 π(r)dr < ∞ を満たすとする．ま
た π(r) が正則変動関数であり，
∫
∞
1
dr
rd−1 π(r)
=∞
(2.14)
を満たすならば，π(kθk) に関する（一般化）ベイズ推定量 δπ は許容的である．
証明
Appendix A．
事前分布の密度関数 π(kθk) に対して，π(r) がインデックス −α の正則変動関数である
とする．このとき，1/{rd−1 π(r)} はインデックス α − d + 1 の正則変動性を持つ．従って，
α − d + 1 > −1 ⇔ α > d − 2
のときに定理 2.3 の条件 (2.14) を満たす．ただし，α > d のときは，狭義の事前分布であ
ることから許容性が従う．一方，
α − d + 1 < −1 ⇔ α < d − 2
のときには (2.14) を満たさない．α = d − 2 が境界であり，(2.14) を満たす関数もあれば，
満たさない関数もある．π(r) = r−(d−2) ，π(r) = r−(d−2) log(r + 1) は (2.14) を満たすが，
どんなに小さい > 0 に対しても π(r) = r−(d−2) {log(r + 1)}1+ は (2.14) を満たさない．
ところで，本稿での事前分布に対する仮定（球面対称性，密度の存在，微分可能性，正
則変動性）のもとで，球面対称な周辺尤度
∫
m̄π (kxk) =
Rd
(
)
1
kx − θk2
exp
−
π(kθk)dθ
2
(2π)d/2
は，Appendix A の補題 A.1 より limr→∞ π(r)/m̄π (r) → 1 を満たす．従って，定理 2.3 は
条件 A2? で条件 A2 を置き換えた定理 2.2 と本質的に同等である．
149
スタイン・パラドクスとベイズ理論
2.2
x の改良
事前分布 ν(dθ) に対する（一般化）ベイズ推定量 x + ∇ log mν (x) が x を改良するため
の十分条件は，(2.5) より
∇T ∇mν (x) k∇mν (x)k2
−
mν (x)
{mν (x)}2
√
∇T ∇ mν (x)
=4 √
≤0
mν (x)
ξ[∇ log mν (x)] = 2
(2.15)
で与えられる．2 回連続微分可能な関数についてラプラシアン ∇T ∇ が任意の x について非
√
正であることは，優調和性として知られている．つまり，周辺尤度の平方根 mν (x) が優
√
調和であれば，x を改良する．また，(2.15) より周辺尤度 mν (x) の優調和性から mν (x)
の優調和性が従う．
さらに，事前分布に対して直接，x を改良するための十分条件を与えよう．事前分布 ν(dθ)
が 2 回連続微分可能な密度関数 π(θ) を持つと仮定する2) ．このとき (2.12) より
}
∫
∫ { 2
exp(−kx − θk2 /2)
∂
exp(−kx − θk2 /2)
∂2
π(θ)dθ
=
π(θ)
dθ
∂x2i Rd
∂θi2
(2π)d/2
(2π)d/2
Rd
{∫
}
∫
exp(−kx − θk2 /2)
exp(−kx − θk2 /2)
T
∇T ∇
π(θ)dθ
=
{∇
∇π(θ)}
dθ
(2π)d/2
(2π)d/2
Rd
Rd
であり，π(θ) の優調和性が mπ (x) の優調和性を含意することが分かる．従って次の定理
が成り立つ．
定理 2.4 (Stein (1981)) 優調和性を持つ事前密度関数 π(θ) に対する（一般化）ベイ
ズ推定量は，x を改良する．
特に，原点に関して球面対称な縮小型事前分布
π(θ) = kθk−α (0 < α < d − 1)
(2.16)
を考え，x を改良するための α に関する十分条件を確認する3) ．kθk−α のラプラシアンは，
θ 6= 0 に対して
∇T ∇π(θ) = α{α − (d − 2)}kθk−α−2
で与えられる．従って，kθk−α は 0 ≤ α ≤ d − 2 のとき優調和性を持ち，従って x を改良
する一般化ベイズ推定量を構成する．ただし 0 < α < d − 2 の場合は，前節の定理 2.3 で与
えられた許容的であるための十分条件を満たさないことが問題として残る．実は，Brown
2)
3)
2.1 節では，密度関数 π は kθk の関数であった．本節では密度関数に球面対称性を仮定せず，π を θ の関数
とする．
α ≥ d − 1 のとき，周辺尤度が原点 x = 0 で発散するため，一般化ベイズ推定量が定義出来ない．
150
日本統計学会誌第45巻第1号 2015
(1971) による一般化ベイズ推定量が非許容的であるための十分条件から，0 < α < d − 2
に対応する一般化ベイズ推定量の非許容性が示される．
√
もちろん π(θ) の優調和性は mπ (x) の優調和性の十分条件の一つに過ぎない．kθk−α は
√
d − 2 < α < d − 1 のとき劣調和的 ∇T ∇π(θ) ≥ 0 である．この場合でも ∇T ∇ mπ (x) ≤ 0
を直接確認することにより，x を改良する一般化ベイズ推定量が見つかる．また，定理 2.3
よりその許容性も従う．
定理 2.5 (Maruyama (1998)，丸山・竹村 (2007)) d − 2 ≤ α ≤ d − 2 + 1/5 とす
る．このとき事前分布 kθk−α に対する一般化ベイズ推定量は許容的であり，また x を改良
する．
証明
3.
Appendix B．
予測分布とパラメータ推定
前節冒頭に述べた統計モデル，すなわち yij (i = 1, . . . , n, j = 1, . . . , d) が互いに独
立な正規分布 N (θj , 1) に従う確率変数を想定する．このとき，同じ分布に従う yij (i =
n + 1, . . . , n + m, j = 1, . . . , d) の予測分布を求める問題を考える．前節と同様に十分
性を考慮し，また簡単のため m = 1 とすると，x ∼ Nd (θ, n−1 Id ) が観測されるときの
z ∼ Nd (θ, Id ) の密度関数 φ(z; θ, 1) を x の関数 p̂(z; x) として与える問題に帰着する．本
節では φ(·; µ, σ 2 ) を d 次元多変量正規分布 Nd (µ, σ 2 Id ) の確率密度関数とする．損失関数
として，φ(z; θ, 1) と p̂(z; x) のカルバック・ライブラー距離を採用し，リスク関数
]
[∫
φ(z; θ, 1)
dz
RKL (θ, p̂) = Ex
φ(z; θ, 1) log
p̂(z; x)
(3.1)
∫
∫
φ(z; θ, 1)
= φ (x; θ, 1/n) φ(z; θ, 1) log
dzdx
p̂(z; x)
で p̂(z; x) の良さを測る．リスク関数 RKL (θ, p̂) のもとで事前分布の密度関数 π(θ) に対す
るベイズ予測分布は，φ(z; θ, 1) の事後分布に関する期待値
∫
φ(z; θ, 1)φ (x; θ, 1/n) π(θ)dθ
∫
p̂π (z; x) =
φ (x; θ, 1/n) π(θ)dθ
で与えられる．
最も単純な予測分布は単に x を θ にプラグインした φ(z; x, 1) であるが，これはベイ
ズ解ではない．Aitchison (1975) は広義の一様分布 π(θ) = 1 に対するベイズ予測分布
p̂U = φ(z; x, (n + 1)/n) がプラグイン予測分布 φ(z; x, 1) を全ての d に対して改良すること
を示した．さらに Komaki (2001) は調和型事前分布 π(θ) = kθk−(d−2) に対するベイズ予
測分布が d ≥ 3 のとき，p̂U を改良することを示した．これは，縮小型事前分布がパラメー
151
スタイン・パラドクスとベイズ理論
タ推定と同様に有効であり，パラメータ推定と予測問題の関連を示唆する結果であった．
George et al. (2006) 及び Brown et al. (2008) は，その関連を明確に与えた．
定理 3.1 (Brown et al. (2008)) 一様分布に関するベイズ予測分布 p̂U と事前分布 π(θ)
に関するベイズ予測分布 p̂π のリスクの差 RKL (θ, p̂U ) − RKL (θ, p̂π ) は，推定問題における一
様分布に関するベイズ推定量と π に関するベイズ推定量のリスクの差を用いて
∫
1/n
RKL (θ, p̂U ) − RKL (θ, p̂π ) =
1/(n+1)
(
Ev [kδU − θk2 ] Ev [kδπ − θk2 ]
−
`
`
)
d`
2`
と表現できる．ただし，v ∼ Nd (θ, `Id ) とし，上記の「推定問題」とは尺度について基準
化された二乗損失関数 kδ − θk2 /` のもとでの θ の推定問題である．δU = v と δπ は，それ
ぞれ一様分布と事前分布の密度関数 π(θ) に対するベイズ推定量である．
証明
Appendix C．
定理 3.1 より，優調和性を持つ事前分布の密度関数 π(θ) に対するベイズ予測分布は，p̂U
を改良する．また，劣調和な場合を含む縮小事前分布 π(θ) = kθk−α（0 < α ≤ d − 2 + 1/5）
に対するベイズ予測分布が p̂U を改良することも，前節の定理 2.5 から従う．
4.
線形回帰モデルのパラメータ推定
目的変数 y を d 個の説明変数 x1 , . . . , xd で説明する線形回帰モデル
y = α1n + Xβ + σ
(4.1)
を考える．ここでサンプルサイズは n とし，(4.1) の各項は n 次元ベクトルである．また α
は未知の切片，1n は各成分が 1 の n 次元ベクトル，X = (x1 , . . . , xd ) は n × d 説明変数行
列，β は d 次元回帰係数ベクトルである．また X の各列は 1 ≤ i ≤ d に対して xTi 1n = 0 と
基準化されているとし，n > d + 1 及び x1 , . . . , xd の一次独立性を仮定する．さらに (4.1)
の誤差項において，σ は未知の尺度パラメータ， = (1 , . . . , n )T は原点周りの球面対称分
布である．その確率密度関数は
∼ f (kk2 )
(4.2)
で与えられ，E[] = 0n 及び Var[] = In を仮定する．従って，誤差項 σ の各成分の期待
値，分散，共分散は，E[σi ] = 0，var(σi ) = σ 2 ，i 6= j について cov(σi , σj ) = 0 を満た
す．なお f (kk2 ) が確率密度関数であること及び Var[] = In は，f に対する条件として
∫
f (kk2 )d =
Rn
π n/2
Γ(n/2)
∫
0
∞
sn/2−1 f (s)ds = 1
(4.3)
152
日本統計学会誌第45巻第1号 2015
および
∫
Rn
kk2 f (kk2 )d =
π n/2
Γ(n/2)
∫
∞
sn/2 f (s)ds = n
(4.4)
0
に対応する．
線形回帰モデルにおいては，多変量正規分布に従う誤差項が通常仮定されるが，本稿で
は多変量 t 分布のように裾が重い分布を含め，f の関数形に出来るだけ条件を置かずに，提
案する推定量の良さを議論する．特に Maruyama (2003)，Maruyama and Strawderman
(2005, 2013) の結果を概説する．
4.1
回帰係数ベクトル β の推定
回帰係数ベクトル β の推定問題を，尺度調整した予測二乗損失関数
L(β̂, β, σ 2 ) =
kX β̂ − Xβk2
σ2
(4.5)
のもとで考える．予測二乗損失関数と呼ばれるのは以下のように説明される．モデル (4.1)
から
y∗ = α1n + Xβ + σ∗
のように新たな y∗ を得るという状況を考える．このとき，期待値 E[y∗ ] = α1n + Xβ を
α̂1n + X β̂ で推定する問題において，尺度調整した二乗損失関数は
n(α̂ − α)2 + kX β̂ − Xβk2
k(α̂1n + X β̂) − (α1n + Xβ)k2
=
σ2
σ2
で与えられる．(4.5) は回帰係数 β に関する項に対応するため予測二乗損失関数と呼ばれ
る．以下では，記述の簡単化のため η = 1/σ 2 ，θ = (X T X)1/2 β 及び θ̂ = (X T X)1/2 β̂ に
対する損失関数
L(θ̂, θ, η) = ηkθ̂ − θk2
(4.6)
のもとでの θ の推定問題として議論を進める．
β の最小二乗推定量
β̂OLS = (X T X)−1 X T y
(4.7)
に対応する
u = (X T X)1/2 β̂OLS = (X T X)−1/2 X T y
は θ の最も標準的な推定量であり，無情報事前分布 π(α, θ, η) = η −1 のもとでの一般化ベ
イズ推定量でもある．しかし第 2 節におけるスタイン・パラドクスと同様に d ≥ 3 のとき，
スタイン・パラドクスとベイズ理論
153
u を改良する推定量が存在する．以下では，f に関して一様に u を改良する一般化ベイズ
推定量を構成していく．事前分布として
π(α, θ, η) = π(θ, η) = kθk−a η b−1
(4.8)
という特別な形をした広義の事前分布を考える．θ に関しては，第 2 節のように縮小型事
前分布を想定するため，a > 0 である．η に関しては，無情報事前分布からのずれを b ∈ R
で調整する．
損失関数 (4.6) のもとでの事前分布 (4.8) に対する θ の一般化ベイズ推定量は，ηkθ̂ − θk2
の事後分布に関する期待値を最小にする E[ηθ | y]/E[η | y] であり，
)
∫∫∫ n/2 (
θη f ηky − α1n − X(X T X)−1/2 θk2 kθk−a η b dαdθdη
E[ηθ | y]
)
= ∫∫∫ n/2 (
E[η | y]
η f ηky − α1n − X(X T X)−1/2 θk2 kθk−a η b dαdθdη
(4.9)
と書ける．ここで (4.9) の分母分子における η に関する積分
∫
(
)
η n/2+b f ηky − α1n − X(X T X)−1/2 θk2 dη
に注目する．
∫∞
0
sn/2+b f (s) ds
<
∞ のとき，正規分布に対応する fG (t)
=
(2π)−n/2 exp(−t/2) に対して，
∫ ∞
(
)
η n/2+b f ηky − α1n − X(X T X)−1/2 θk2 dη
0
∫ ∞
= ky − α1n − X(X T X)−1/2 θk−2(n/2+b+1)
sn/2+b f (s) ds
(4.10)
0
∫ ∞ n/2+b
∫ ∞
(
)
s
f (s) ds
η n/2+b fG ηky − α1n − X(X T X)−1/2 θk2 dη
= ∫ 0∞ n/2+b
s
fG (s) ds 0
0
∫∞
∫∞
を得る．Maruyama (2003) は定数 0 sn/2+b f (s) ds/ 0 sn/2+b fG (s) ds が分母分子でキャ
ンセルされるので，事前分布 (4.8) のもとでの一般化ベイズ推定量 (4.9) は，f の関数形に
は依存せず，正規分布のもとでの一般化ベイズ推定量
(
)
∫∫∫ n/2+b
θη
exp −ηky − α1n − X(X T X)−1/2 θk2 /2 kθk−a dαdθdη
E[ηθ | y]
(
)
= ∫∫∫ n/2+b
E[η | y]
η
exp −ηky − α1n − X(X T X)−1/2 θk2 /2 kθk−a dαdθdη
(4.11)
∫ ∞ n/2+b
に一致することを示した．なお 0 s
f (s) ds < ∞ は誤差項 i の 2 + 2b 次のモーメン
トの存在に対応する．従って，多変量正規分布のもとでの一般化ベイズ推定量 (4.11) が一
般の球面対称分布でも一般化ベイズ推定量であることの正当化には，b > 0 の場合に 2 次
より大きい (2 + 2b) 次モーメントの存在が必要となる．
さて，ky − α1n − X(X T X)−1/2 θk2 を，最小二乗推定量 u と残差平方和
z = k(In − X(X T X)−1 X T )(y − ȳ1n )k2
154
日本統計学会誌第45巻第1号 2015
によって
ky − α1n − X(X T X)−1/2 θk2 = n(α − ȳ)2 + ku − θk2 + z
(4.12)
と分解する．すると (4.11) と (4.12) より，事前分布 (4.8) に関する周辺尤度
(
)
∫∫∫
2
2
n/2+b−1 exp −η{n(α − ȳ) + ku − θk + z}/2
mπ (u, z) =
η
kθk−a dαdθdη
(2π)n/2
に対して，一般化ベイズ推定量は
1 ∇u mπ (u, z)
E[ηθ | y]
=u−
E[η | y]
2 {∂/∂z}mπ (u, z)
(4.13)
と表せる4) ．mπ (u, z) の α，θ ，η に関する積分を行うにあたり，u と z に依存しない定数
は (4.13) の分母分子でキャンセルされるため，以下の積分計算では単に c1 , c2 , c3 , c4 のよう
に記す．α に関する積分により
∫∫
(
)
mπ (u, z) = c1
η (n−1)/2+b−1 exp −η{ku − θk2 + z}/2 kθk−a dθdη
を得る．また θ に関する積分において，
∫∞
0
ta/2−1 exp(−tηkθk2 /2)dt ∝ η −a/2 kθk−a 及び
平方完成
tkθk2 + kθ − uk2 = (1 + t)kθ − u/(1 + t)k2 + tkuk2 /(1 + t)
に注意すると
∫
∫
(
)
exp −ηku − θk2 /2 kθk−a dθ = c2 η (a−d)/2
0
∞
(
)
ta/2−1
ηtkuk2
exp
−
dt
2(1 + t)
(1 + t)d/2
を得る．さらに η に関する積分は
(n − 1 + a − d)/2 + b > 0
(4.14)
のとき
∫
∞
(
tkuk2 + z(1 + t)
η
exp −η
2(1 + t)
0
(
)
−(n−1+a−d)/2−b
t{kuk2 + z} + z
= c3
1+t
(n−1+a−d)/2+b−1
)
dη
である．まとめると，mπ (u, z) は t による積分で表現されて
∫
∞
mπ (u, z) = c4
0
4)
ta/2−1 (1 + t)(n−1+a−d)/2+b−d/2
dt
(t{kuk2 + z} + z)(n−1+a−d)/2+b
(4.15)
この表現は，第 2 節で (2.7) によって与えられたベイズ推定量の周辺尤度による表現 x + ∇ log mν (x) に対
応している．
155
スタイン・パラドクスとベイズ理論
である．(4.13) に (4.15) を適用すると，θ の一般化ベイズ推定量
(
)
∫ ∞ a/2
t (1 + t)(n−1+a−d)/2+b−d/2 (t{kuk2 + z} + z)−(n−1+a−d)/2−b−1 dt
0
1 − ∫ ∞ a/2−1
u
t
(1 + t)(n−1+a−d)/2+b−d/2+1 (t{kuk2 + z} + z)−(n−1+a−d)/2−b−1 dt
0
(4.16)
を得る．
Maruyama and Strawderman (2005) は，
(n − 1 + a − d)/2 + b − d/2 = 0
(4.17)
の場合に，(4.15) において t に関する積分が実行でき，mπ (u, z) 及び一般化ベイズ推定量
が簡潔な表現を持つことを示した．実際
∫ ∞
mπ (u, z) = c4
ta/2−1 (t{kuk2 + z} + z)−d/2 dt
0
∫ ∞
= c4 z −d/2
ta/2−1 (t{kuk2 /z + 1} + 1)−d/2 dt
0
∫ ∞
−d/2
= c4 z
(kuk2 /z + 1)−a/2
sa/2−1 (s + 1)−d/2 ds
0
= c5 z −d/2 (kuk2 /z + 1)−a/2
と計算される．また ν = a/(d − a) を導入すると，一般化ベイズ推定量
(
)
ν
1 ∇u mπ (u, z)
= 1−
u
θ̂S = u −
2 {∂/∂z}mπ (u, z)
kuk2 /z + 1 + ν
(4.18)
を得る．なお，a ∈ (0, d) より ν ∈ (0, ∞) であり，(4.17) は ν を用いると
b = − {(n − 1)/2 − d} − (d/2){ν/(ν + 1)}
(4.19)
と書ける．従って，d ≤ (n − 1)/2 ならば，任意の ν ∈ (0, ∞) に対して (4.17) を満たす b
は b < 0 を満たす．つまり，Var[] = In のもとで (4.18) がベイズ推定量として正当化出来
る．一方 d > (n − 1)/2 の場合，ν ≤ (2d − n + 1)/(n − d − 1) なる小さい ν について，θ̂S
のベイズ性を正当化するには，2 次より大きい高次モーメントの存在が必要である．
さらに，θ̂S が最小二乗推定量 u を改良するための十分条件について考える．縮小推定量
(1 − ψ(kuk2 /z)/{kuk2 /z})u が u を改良するための十分条件として，が 0n に関して球
面対称な誤差項であり Var[] = In という仮定のもとで，
ψ が単調非減少で絶対連続， 0 ≤ ψ ≤ 2(d − 2)/(n − d + 1)
が知られている5) ．従って 0 < ν ≤ 2(d − 2)/(n − d + 1) のとき，θ̂S は u を改良する．
5)
証明については Maruyama and Strawderman (2005) の Theorem 2.1 を参照されたい．
156
日本統計学会誌第45巻第1号 2015
以上の議論を，もとのパラメータ β に戻し，また kuk2 /z + 1 が回帰分析の決定係数 r2
を用いて 1/(1 − r2 ) と書けることに注意して，θ̂S に対応する推定量の性質を定理の形にま
とめる．
定理 4.1 (Maruyama and Strawderman (2005)) 線形回帰モデル
y = α1n + Xβ + σ
において，β の推定問題を損失関数 kX β̂ − Xβk2 /σ 2 のもとで考える．は 0n に関して球
面対称な誤差項であり，Var[] = In を仮定する．このとき，ν ∈ (0, ∞) 及び決定係数 r2
を用いた，最小二乗推定量 β̂OLS に対する縮小推定量
β̂S =
β̂OLS
ν(1 − r2 ) + 1
は以下の性質を持つ．
1. 0 < ν ≤ 2(d − 2)/(n − d + 1) のとき，β̂OLS を改良する．
2. max(0, (2d − n + 1)/(n − d − 1)) < ν < ∞ に対して，一般化ベイズ推定量である．
3. d > (n − 1)/2 のとき 0 < ν < (2d − n + 1)/(n − d − 1) に対して，の {2 + (2d −
n + 1) − dν/(ν + 1)} 次モーメントの存在のもとで一般化ベイズ推定量である．
Remark 4.1
(
)
θ̂S = 1 − ν/{kuk2 /z + 1 + ν} u は，スタイン・パラドクスが見出さ
れた Stein (1956) で考えられた推定量に関係している．Stein (1956) は，第 2 節における
x ∼ Nd (θ, Id ) に対する θ の推定問題で δa,b = (1 − a/{a + b + kxk2 })x という推定量を考
えて，小さい a > 0 と大きい b > 0 に対して，x を改良することを示した．θ̂S は δa,b の分
散未知の場合に対応することが分かる．
その後 Strawderman (1971) をはじめ多くの研究者によって，第 2 節における x，本節
における u を改良する（一般化）ベイズ推定量が提案されてきた．しかし，それらは全て
縮小係数が (4.16) のように積分の比で表されており，x や u を改良するベイズ推定量はそ
のように複雑なものだと考えられてきた．簡潔な形をもつ推定量 θ̂S の一般化ベイズ推定量
としての正当化はその意味で興味深い．
ところで，Stein (1956) の推定量 δa,b = (1 − a/{a + b + kxk2 })x は，
ma,b (kxk) = (kxk2 + a + b)−a/2
に対して，x + ∇ log ma,b (kxk) と表現できる．従って，
∫
1
(kxk2 + a + b)−a/2 =
exp(−kx − θk2 /2)ν(dθ)
d/2
Rd (2π)
(4.20)
スタイン・パラドクスとベイズ理論
157
なる非負の測度 ν があれば，δa,b は一般化ベイズ推定量として解釈出来る．しかし，Maruyama
and Iwasaki (2005) の方法を用いると，(4.20) を満たす ν は符号付き測度になるはずであ
る．従って，定理 2.1 より δa,b は非許容的であると予想される．
4.2
分散 σ 2 の推定
誤差項の分散 σ 2 = E[{σi }2 ] の推定問題をスタイン損失関数
LS (δ, σ 2 ) = δ/σ 2 − log(δ/σ 2 ) − 1
に関するリスク R({α, θ, σ 2 }, δ) = E[LS (δ, σ 2 )] の下で考える．正規分布の仮定のもとで，
不偏推定量
δU = z/(n − d − 1)
の非許容性が Stein (1964) によって示された．具体的には，β = θ = 0 の場合の σ 2 の不
偏推定量 ky − ȳ1n k2 /(n − 1) と δU との最小値をとる推定量
(
)
ky − ȳ1n k2
δST = min δU ,
n−1
が δU を改良することが示された．Brewster and Zidek (1974) は ψ(r2 )δU が δU を改良する
ための十分条件
1. ψ(·) が単調非減少
2. ψBZ (r2 ) ≤ ψ(r2 ) ≤ 1，ただし ψBZ (r2 ) は滑らかな増加関数であり，
ψBZ (r2 ) = 1 −
2(1 − r2 )(n−d−1)/2
∫1
(n − 1) 0 td/2−1 (1 − r2 t)(n−d−1)/2 dt
を導出した．ψBZ 自身および ψST はこの十分条件を満たす．また，Ghosh (1994) は Brewster
and Zidek (1974) の十分条件を満たす一般化ベイズ推定量のクラスを与えた．そのクラス
の中には kθk−(d−2) η −1 に対する一般化ベイズ推定量が含まれる．
以下では，球面対称分布に従う誤差項の分散推定の結果を紹介する．Maruyama and
Strawderman (2013) は，事前分布
π(α, θ, η) = π(α, θ)η −1
に対する一般化ベイズ推定量 1/E[η | y] が，f に関して一様に，正規分布のもとでの一
般化ベイズ推定量に一致することを示した．これは次のように確認出来る．mi (y | α, θ)
(i = 0, 1) を y の α 及び θ を与えた時の η i−1 に関する条件付き周辺尤度
∫ ∞
(
)
mi (y | α, θ) =
η n/2 f η{n(α − ȳ)2 + kθ̂ − θk2 + z} η i−1 dη
0
158
日本統計学会誌第45巻第1号 2015
とする．このとき一般化ベイズ推定量は
∫∫
m (y | α, θ)π(α, θ)dαdθ
∫∫ 0
m1 (y | α, θ)π(α, θ)dαdθ
と書ける．また (4.10) と同様に fG (s) = (2π)−n/2 exp(−s/2) に対して
∫ ∞
mi (y | α, θ) = {n(α − ȳ)2 + kθ̂ − θk2 + z}−n/2−i
s{n+2i}/2−1 f (s)ds
0
∫ ∞ (n+2i)/2−1
∫
)
s
f (s)ds ∞ n/2 (
= ∫ 0∞ (n+2i)/2−1
η fG η{n(α − ȳ)2 + kθ̂ − θk2 + z} η i−1 dη
s
fG (s)ds 0
∫0∞ (n+2i)/2−1
s
f (s)ds G
= ∫ 0∞ (n+2i)/2−1
mi (y | α, θ)
s
fG (s)ds
0
である．従って，一般化ベイズ推定量は
∫∫
m (y | α, θ)π(α, θ)dαdθ
∫∫ 0
m1 (y | α, θ)π(α, θ)dαdθ
∫∞
∫ ∞ n/2−1
∫∫ G
s
f (s)ds 0 sn/2 fG (s)ds
m (y | α, θ)π(α, θ)dαdθ
0∫
∫∫ G0
∫∞
=
∞ n/2
n/2−1
m1 (y | α, θ)π(α, θ)dαdθ
s f (s)ds 0 s
fG (s)ds
∫∫0 G
m (y | α, θ)π(α, θ)dαdθ
= ∫∫ 0G
m1 (y | α, θ)π(α, θ)dαdθ
となり，正規分布のもとでの一般化ベイズ推定量に一致する．ただし，最後の等式で (4.3)
と (4.4) を用いた．
さらに Maruyama and Strawderman (2013) は，正規分布のもとで δU を改良する ψ(r2 )δU
は，球面対称分布の部分クラスである正規分布の尺度混合分布
∫
∞
f (t) =
(2πτ 2 )−n/2 exp(−t/{2τ 2 })g(τ 2 )dτ 2
0
（ただし τ は E[τ 2 ] = 1 を満たす）のもとでも δU を改良することを示した．
以上をまとめると，事前分布 π(α, θ, η) = kθk−(d−2) η −1 に対する一般化ベイズ推定量の
良さが以下のように主張できる．
定理 4.2 (Maruyama and Strawderman (2013))
事前分布 π(α, θ, η)
=
kθk−(d−2) η −1 を想定する．このとき，球面対称分布のもとでの一般化ベイズ推定量は，
正規分布のもとでの一般化ベイズ推定量に一致する．また，その一般化ベイズ推定量は，
正規分布の尺度混合分布に対して一様に，不偏推定量 δU を改良する．
事前分布 π(α, θ, η) = kθk−(d−2) η −1 に対する一般化ベイズ推定量の具体的な形などの詳
細は Maruyama and Strawderman (2013) を参照されたい．
スタイン・パラドクスとベイズ理論
Appendix A.
A.1
159
定理 2.3 の証明
Blyth の十分条件
一般化ベイズ推定量が許容的であるための十分条件として，Blyth (1951) の方法が知ら
れている．いくつかのバージョンがあるが，我々の目的には，Brown (1971) によるものが
使いやすい．
定理 A.1 (Blyth (1951)，Brown (1971)) 広義の事前分布の密度関数 π(θ) が
∫
∫
π(θ)dθ = ∞,
π(θ)dθ = cπ > 0
Rd
kθk≤1
を満たすとする．また任意に固定した θ に対して，π に収束する増加列 0 ≤ π1 (θ) ≤ π2 (θ) ≤
· · · ≤ π(θ) は，kθk ≤ 1 の領域で任意の i について πi (θ) ≡ π(θ) を満たし，また任意に固
∫
定した i について Rd πi (θ)dθ < ∞ を満たすとする（ただし必ずしも積分値が 1 にはなら
ない）．このとき，πi (θ) に対するベイズ推定量 δπi と π(θ) に対する一般化ベイズ推定量
δπ の πi (θ) に対する（基準化されない）ベイズリスクの差
∫
∆i = {R(θ, δπ ) − R(θ, δπi )} πi (θ)dθ
が limi→∞ ∆i = 0 を満たすならば，一般化ベイズ推定量 δπ は許容的である．
証明
背理法の仮定として，δπ が非許容的であり，R(θ, δ 0 ) ≤ R(θ, δπ ) を満たす δ 0 が
存在するとする．このとき推定量 δ∗ = (δπ + δ 0 )/2 は Jensen の不等式より，任意の θ に
ついて
R(θ, δ∗ ) < [R(θ, δ 0 ) + R(θ, δπ )] /2 ≤ R(θ, δπ )
である．ここで R(θ, δ∗ ) と R(θ, δπ ) はその構成方法から θ について連続であり，kθk ≤ 1
に対して R(θ, δ∗ ) < R(θ, δπ ) − を満たすような > 0 が存在する．従って，
∫
∆i ≥
[R(θ, δπ ) − R(θ, δ∗ )]πi (θ)dθ
d
∫R
≥
[R(θ, δπ ) − R(θ, δ∗ )]π(θ)dθ
kθk≤1
≥ cπ > 0
である．これは ∆i → 0 に矛盾する．
A.2
π への収束列
以下では，広義の事前分布の密度関数が，球面対称で微分可能性を持つ π(kθk) である場
合に，定理 A.1 における関数列 π1 , π2 , . . . を構成する．定理 2.3 の仮定 (2.14) すなわち
∫ ∞
du
=∞
d−1
u π(u)
1
160
日本統計学会誌第45巻第1号 2015
のもとで hi (r) を




1


∫r


1/{ud−1 π(u)}du
1
hi (r) = 1 − ∫ i+1

1/{ud−1 π(u)}du

1




0
0≤r≤1
1<r <i+1
(A.1)
r ≥i+1
と定義する．hi (r) は固定した r に対し i の増加関数であり，また limi→∞ hi (r) = 1 であ
る．従って
πi (r) = π(r)h2i (r)
(A.2)
は，固定した r に対し i に関する増加列 π1 (r) ≤ π2 (r) ≤ · · · であり，π(r) に収束する．ま
た 0 ≤ r ≤ 1 において πi (r) = π(r) である．さらに r = i + 1 で打ち切られているから
∫
π (kθk)dθ < ∞ を満たす．
Rd i
hi (r) は，r に関して区分的に微分可能であり，



0




1/{rd−1 π(r)}
d
hi (r) = − ∫ i+1
d−1

dr

 1 1/{u π(u)}du



0
0≤r≤1
1<r <i+1
r ≥i+1
となる．r = 1 及び r = i + 1 では微分不可能であるが，後の議論には影響しない．また固
定した r について，
lim {d/dr}hi (r) = 0
(A.3)
i→∞
が従う．さらに ({d/dr}hi (r))2 を i に依存せず上から抑える関数として



0



)2 
(
 1/{r2(d−1) π 2 (r)}
d
∫2
hi (r) =
sup
( 1 1/{ud−1 π(u)}du)2

dr
i



1/{r2(d−1) π 2 (r)}


 ∫r
( 1 1/{ud−1 π(u)}du)2
を得る．
0≤r≤1
1<r<2
r≥2
(A.4)
スタイン・パラドクスとベイズ理論
A.3
161
ベイズリスクの差
δπi を狭義の事前分布の密度関数 πi (kθk) = π(kθk)h2i (kθk) に対するベイズ推定量とす
る．このとき，δπ と δπi の πi (kθk) に関する基準化されないベイズリスクの差は
∫
[R(θ, δπ ) − R(θ, δπi )] πi (kθk)dθ
∫ ∫
=
[kδπ − θk2 − kδπi − θk2 ]φ(x; θ, 1)πi (kθk)dθdx
Rd Rd
}
∫ {
∫
2
2
T
=
[kδπ k − kδπi k ]m(πi ; x) − 2(δπ − δπi )
θφ(x; θ, 1)πi (kθk)dθ dx
d
Rd
∫R
kδπ − δπi k2 m(πi ; x)dx
=
∆i =
Rd
Rd
と表現できる．ただし φ(x; µ, σ 2 ) は d 次元多変量正規分布 Nd (µ, σ 2 Id ) の確率密度関数で
あり，また m(g(θ); x) は θ ∼ Nd (x, Id ) に対する g(θ) の期待値である．δπi の構成方法か
ら，各点 x で limi→∞ δπi (x) = δπ (x) であり，kδπ − δπi k2 m(πi ; x) を上から可積分関数で
抑えれば優収束定理より ∆i → 0 が従う．
ここでベイズ推定量の表現 (2.11)，(2.12) を用いると
m(∇π; x) m(∇{πh2i }; x) 2
2
∆i =
m(π; x) − m(πh2 ; x) m(πhi ; x)dx
Rd
i
∫ m(∇π; x) m({∇π}h2i ; x) m(π∇h2i ; x) 2
2
−
−
=
2 ; x)
2 ; x) m(πhi ; x)dx
m(π; x)
m(πh
m(πh
d
R
i
i
∫
を得る．さらに三角不等式より ∆i ≤ 2(∆1i + ∆2i ) を得る．ここで
m(π∇h2i ; x) 2
2
∆1i =
m(πh2 ; x) m(πhi ; x)dx
Rd
i
∫ m(∇π; x) m({∇π}h2i ; x) 2
m(πh2i ; x)dx
∆2i =
−
2 ; x)
m(π; x)
m(πh
d
R
i
∫
である．A.4 節と A.5 節で，それぞれ ∆1i と ∆2i を可積分関数で上から抑える．
A.4
∆1i
コーシー・シュワルツの不等式より
∫
∆1i = 4
2
∫R
d
≤4
∫
Rd
=4
Rd
km(πhi ∇hi ; x)k
dx
m(πh2i ; x)
m(πk∇hi k2 ; x)dx
πk∇hi k2 dθ = 4
2π d/2
Γ(d/2)
∫
∞
0
rd−1 π(r){h0i (r)}2 dr
162
日本統計学会誌第45巻第1号 2015
を得る．また (A.4) より
∫ ∞
Γ(d/2)
∆1i ≤
rd−1 π(r) sup{h0i (r)}2 dr 8π d/2
i
0
∫ 2
∫ ∞
1/{rd−1 π(r)}
1/{rd−1 π(r)}
∫r
=
dr +
∫2
( 1 1/{ud−1 π(u)}du)2
1 ( 1 1/{ud−1 π(u)}du)2
2
[
]∞
1
1
+ −∫ r
= ∫2
1/{ud−1 π(u)}du
1/{ud−1 π(u)}du
1
1
= ∫2
1
2
2
1/{ud−1 π(u)}du
<∞
である．
A.5
∆2i
∆2i の被積分関数の θj に関する偏微分について
(
)2
m(∇j π; x) m({∇j π}h2i ; x)
−
m(π; x)
m(πh2i ; x)
(
)2
1
m(∇j π; x)
2
2
m(πhi ; x)
− m({∇j π}hi ; x)
m(πh2i ; x)
m(π; x)
( (
{
} ))2
1
m(∇j π; x) ∇j π
2
m
πh
−
;x
i
m(πh2i ; x)
m(π; x)
π
(
}2 )
{
m(πh2i /k; x)
m(∇j π; x) ∇j π
2
m πhi k
−
;x
m(πh2i ; x)
m(π; x)
π
( {
}2 )
m(πh2i /k; x)
m(∇j π; x) ∇j π
m πk
−
;x
m(πh2i ; x)
m(π; x)
π
m(πh2i ; x)
=
=
≤
≤
(A.5)
を得る．ただし
k(kθk) =


kθk
0 < kθk < 1

1
otherwise
である．(A.5) における最初の不等式はコーシー・シュワルツの不等式から，また 2 番目の
不等式は h2i ≤ 1 から従う．
(A.5) において m(πh2i /k; x)/m(πh2i ; x) は，
m({π/kθk}I[0,1] (kθk); x) + m(πh2i I(1,∞) (kθk); x)
m(πh2i /k; x)
=
2
m(πhi ; x)
m(πI[0,1] (kθk); x) + m(πh2i I(1,∞) (kθk); x)
m({π/kθk}I[0,1] (kθk); x)
≤
m(πI[0,1] (kθk); x)
(A.6)
と評価できる．補題 A.2 から，(A.6) の右辺は i および x に依存しない定数で上から抑え
られる．
163
スタイン・パラドクスとベイズ理論
ついで (A.5) における
(
{
m πk
m(∇j π; x) ∇j π
−
m(π; x)
π
)
}2
;x
(A.7)
は
)2
(
)
k{∇j π}2
m(∇j π; x)
m(∇j π; x)
+m
; x − 2m (k∇j π; x)
m(πk; x)
m(π; x)
π
m(π; x)
(
)
m(∇j π; x)2
k{∇j π}2
m(∇j π; x)
≤
+m
; x − 2m (k∇j π; x)
m(π; x)
π
m(π; x)
)
(
m(∇j π; x)2
k{∇j π}2
m(∇j π; x)
;x −
+ 2m ({1 − k}∇j π; x)
=m
π
m(π; x)
m(π; x)
(
(A.8)
と評価できる．(A.8) の右辺の第 3 項に関しては，補題 A.3 より可積分である．また，j =
1, . . . , d で右辺の第 3 項に関する和をとっても可積分である．さらに，(A.8) の最右辺の第
1 項，第 2 項を j = 1, . . . , d に関して和を取った
(
m
) ∑
d
kk∇πk2
m(∇j π; x)2
;x −
π
m(π; x)
j=1
(A.9)
は，補題 A.4 より可積分である．
A.6
補題
補題 A.1 (Maruyama and Takemura (2008)) y = (y1 , . . . , yd )T ∼ Nd (µ, Id ) と
する．%(r) は微分可能であり，r ≥ r∗ のとき %(r) が符号を変えないような r∗ が存在する
ことを仮定する．さらに正則変動性を持つとする．このとき
E[%(kyk)]
E[yi %(kyk)]
− 1 及び − µi kµk %(kµk)
%(kµk)
は kµk ≥ r∗ で有界である．
補題 A.2
∫1
m({π/kθk}I[0,1] (kθk); x)
π(r)rd−2 dr
< exp(1/2) ∫01
m(πI[0,1] (kθk); x)
π(r)rd−1 dr
0
証明
kθk2 が自由度 d，非心度パラメータ kxk2 の非心カイ二乗分布に従うので，
aj (x) =
1
(kxk2 /2)j
exp(−kxk2 /2)
d/2+j
j!
Γ(d/2 + j)2
に対して，
m({π/kθk}I[0,1] (kθk); x)
=
m(πI[0,1] (kθk); x)
=
∫1
d−1+2j
exp(−r2 /2)dr
j=0 aj (x) 0 {π(r)/r}r
∫1
∑∞
d−1+2j exp(−r 2 /2)dr
j=0 aj (x) 0 π(r)r
∑∞
aj (x)E[r2j−1 ]
∑j=0
∞
2j
j=0 aj (x)E[r ]
∑∞
164
日本統計学会誌第45巻第1号 2015
を得る．ここで最右辺は，確率密度関数
π(r)rd−1 exp(−r2 /2)I[0,1] (r)
∫1
π(r)rd−1 exp(−r2 /2)dr
0
に従う確率変数 r に対する期待値である．期待値の比に関して，
E[r−1 ] ≥
より，
E[r]
E[r3 ]
≥
≥ ···
E[r2 ]
E[r4 ]
∑∞
aj (x)E[r2j−1 ]
∑j=0
∞
2j
j=0 aj (x)E[r ]
∫1
≤ E[r
−1
] = ∫01
0
π(r)rd−2 exp(−r2 /2)dr
π(r)rd−1 exp(−r2 /2)dr
を得る．0 ≤ r ≤ 1 に対して，exp(−1/2) ≤ exp(−r2 /2) ≤ 1 であるから，補題を得る．
∫
補題 A.3
m ({1 − k}∇j π; x) m(∇j π; x) dx < ∞
m(π; x) d
R
証明
まず |∇j π| ≤ |π 0 | に注意する．また定理の仮定より rπ 0 (r)/π(r) が原点付近で有
界であり，さらに π(r) は正則変動性を有するので rπ 0 (r)/π(r) は r → ∞ で r のインデッ
クスに収束する．よって，r > 0 において |rπ 0 (r)/π(r)| は有界であり，
|rπ 0 (r)/π(r)| < C0
なる C0 が存在する．従って，被積分関数について
(
)
m ({1 − k}∇j π; x) m(∇j π; x) ≤ C02 m {π/kθk}I[0,1] (kθk); x m(π/kθk; x)
m(π; x)
m(π; x)
(A.10)
(A.11)
を得る．以下では，領域 kxk ≤ 1 と kxk ≥ 1 のそれぞれについて (A.11) の右辺の可積分
∫1
性を検討する．なお定理 2.3 の正則条件 0 rd−2 π(r)dr < ∞ はヤコビアンを考慮すると
∫
kθk≤1
π(kθk)
dθ < ∞
kθk
(A.12)
に対応することに注意する．
領域 kxk ≤ 1 での可積分性について，(kθk − 1)2 ≤ kx − θk2 ≤ (kθk + 1)2 に注意すると
}
{
(
) m(π/kθk; x)
m(π/kθk; x)2
sup m {π/kθk}I[0,1] (kθk); x
≤ sup
m(π; x)
m(π; x)
kxk≤1
kxk≤1
)
(∫
2
{π(kθk)/kθk} exp(−{kθk − 1}2 /2)dθ
1
Rd ∫
≤
π(kθk) exp(−{kθk + 1}2 /2)dθ
(2π)d/2
Rd
(A.13)
を得る．(A.13) の右辺の分母分子について，kθk ≤ 1 での可積分性は (A.12) から，kθk ≥ 1
での可積分性は π(r) の正則変動性から従う．従って，(A.11) の右辺は kxk ≤ 1 において
可積分である．
スタイン・パラドクスとベイズ理論
165
(A.11) の右辺の領域 kxk ≥ 1 での可積分性について，補題 A.1 より kxk → ∞ のとき
kxkm(π/kθk; x)/m(π; x) → 1 である．また，minkθk≤1 kx − θk2 = (kxk − 1)2 より，
∫
(
) exp(−{kxk − 1}2 /2)
π(kθk)
dθ
m {π/kθk}I[0,1] (kθk); x ≤
(2π)d/2
kθk≤1 kθk
を得る．従って (A.11) の右辺は kxk → ∞ のとき指数関数のオーダーで減衰するので，
kxk ≥ 1 での可積分性が示された．

 ∫
(
) ∑
d
2
2
kk∇πk
m(∇
π;
x)
j
m
 dx < ∞
;x −
補題 A.4 π
m(π; x)
Rd
j=1
証明
∫
(A.9) の第 1 項の x に関する積分は，積分の順序交換により，
(
)
∫
∫
kk∇πk2
k(kθk)k∇π(kθk)k2
k(kxk)k∇π(kxk)k2
; x dx =
dθ =
dx
m
π
π(kθk)
π(kxk)
Rd
Rd
Rd
となる．従って，(A.9) の代わりに，被積分関数
d
2
k(kxk)k∇π(kxk)k2 ∑
m(∇j π; x) −
π(kxk)
m(π; x) (A.14)
j=1
に関する可積分性を議論すればよい．
また，定理の仮定を満たす π(r) はインデックスが負の正則変動関数なので，ある r0 が
あって，r ≥ r0 に対して π 0 (r) < 0 である．従って，補題 A.1 より kxk ≥ r0 に対して，
0
m(∇j π; x) − xj π 0 (kxk) < C1 |π (kxk)| < C0 C1 π(kxk)
(A.15)
kxk
kxk
kxk2
π(kxk)
|m(π; x) − π(kxk)| < C2
(A.16)
kxk
を満たす C1 > 0, C2 > 0 が存在する．ただし，(A.15) の最初の不等号には π 0 の微分可能性
つまり π の二回微分可能性の仮定が必要である．また，(A.15) の 2 番目の不等号は (A.10)
から従う．
(A.14) の kxk ≤ r0 での可積分性を見る．第 1 項については，(A.10) より
(
)2
{π 0 (kxk)}2
π(kxk)
π 0 (kxk)
π(kxk)
k(kxk)k∇π(kxk)k2
≤ kxk
≤
sup
kxk
≤ C02
π(kxk)
π(kxk)
kxk r∈(0,r0 )
π(kxk)
kxk
であるから，(A.12) より kxk ≤ r0 の領域での可積分性が従う．第 2 項について，被積分
関数
M (π; x) =
d
∑
m(∇j π; x)2
j=1
m(π; x)
の各項の分子
(
)
(
)
∫
kx − θk2
π 0 (kθk)
kx − θk2
m(∇j π; x) =
∇j π exp −
dθ =
θj
exp −
dθ
2
kθk
2
Rd
Rd
∫
166
日本統計学会誌第45巻第1号 2015
が，x = 0 のとき 0 になる．従って被積分関数 M (π; x) も M (π; 0) = 0 を満たす．M (π; x)
の連続性から有界閉集合 kxk ≤ r0 において最大値 M∗ が存在して，
∫
∫
kxk≤r0
M (π; x)dx ≤ M∗
dx
kxk≤r0
である．従って，(A.14) は kxk ≤ r0 の領域で可積分である．
次に (A.14) の kxk ≥ r0 での可積分性を見る．(A.15) と (A.16) より，kxk ≥ r0 の領域
において，C3 > 0 が存在して
d
2
k(kxk)k∇π(kxk)k2 ∑
m(∇
π;
x)
j
< C3 π(kxk)
−
π(kxk)
m(π;
x)
kxk3
(A.17)
j=1
のように被積分関数 (A.14) を上から抑えることが出来る．変換のヤコビアン rd−1 を考慮
∫∞
すると，(A.17) の右辺の可積分性は積分 r0 rd−4 π(r)dr の可積分性に帰着する．定理の条
件 (2.14) よりインデックスが −(d − 2) 以下の正則変動関数 π(r) を考えているので，正則
変動関数 rd−4 π(r) のインデックスは −2 以下である．従って (A.14) は kxk ≥ r0 において
可積分である．
Appendix B.
定理 2.5 の証明
事前分布の密度関数 π(θ) = kθk−α に対する周辺尤度は
∫
(
)
exp −kθ − xk2 /2 kθk−α dθ
d
{∫ ∞
}
∫R
(
)
∝
exp −kθ − xk2 /2
tα/2−1 exp(−tkθk2 /2)dt dθ
Rd
0
(
)
∫ ∞
tkxk2
∝
tα/2−1 (t + 1)−d/2 exp −
dt
2(t + 1)
0
(
)
∫ ∞
kxk2
= exp(−kxk2 /2)
tα/2−1 (t + 1)−d/2 exp
dt
2(t + 1)
0
mπ (x) =
(B.1)
と表現できる．ここで指数関数 exp(kxk2 /{2(t + 1)}) をマクローリン展開し項別積分する
ことにより，(B.1) の最右辺は，合流型超幾何関数
M (a, b, z) = 1 +
∞
∑
a · (a + 1) · · · (a + n − 1) z n
n=1
b · (b + 1) · · · (b + n − 1) n!
を用いて，mπ (x) ∝ exp(−kxk2 /2)M ({d − α}/2, d/2, kxk2 /2) と表せる．(2.15) で与えら
れた一般化ベイズ推定量が x を改良するための十分条件
2
∇T ∇mπ (x) k∇mπ (x)k2
−
≤0
mπ (x)
{mπ (x)}2
スタイン・パラドクスとベイズ理論
167
は，合流型超幾何関数に関する隣接等式
bM (a, b, z) − bM (a − 1, b, z) − zM (a, b + 1, z) = 0
(B.2)
(1 + a − b)M (a, b, z) − aM (a + 1, b, z) + (b − 1)M (a, b − 1, z) = 0
(B.3)
を用いて，また z = kxk2 /2 として
M ({d − 2 − α}/2, d/2, z) α M ({d − 2 − α}/2, d/2, z)
α
+d
−
≥0
2
M ({d − α}/2, d/2 + 1, z)
2 M ({d − α}/2, d/2, z)
(B.4)
と書き直せる．ここで d − 2 < α < d − 1 に対して，M ({d − 2 − α}/2, d/2, z) の引数
が {d − 2 − α}/2 ∈ (−1, 0) であることに注意する．従って M ({d − 2 − α}/2, d/2, z) は
M ({d − 2 − α}/2, d/2, 0) = 1 から M ({d − 2 − α}/2, d/2, ∞) = −∞ まで単調に減少する．
M ({d − 2 − α}/2, d/2, z) ≥ 0 なる z については，(B.4) の左辺の第 2 項と第 3 項において，
d > α/2 及び M ({d − α}/2, d/2, z) ≥ M ({d − α}/2, d/2 + 1, z) より，(B.4) の左辺は α/2
より大きい．すなわち，常に (B.4) が成り立つ．一方 M ({d − 2 − α}/2, d/2, z) < 0 なる z
について，−1 < a < 0, b > 0, z ≥ 0 に対して成り立つ不等式
(a + 1)bM (a, b, z) − a(b + 1)M (a + 1, b + 1, z) ≥ 0
を用いると，(B.4) の左辺が下から α/2 − (d + 2)(α − d + 2)/(d − α) で抑えられる．従って，
α(d − α) + 2(d + 2)(d − α − 2) ≥ 0
であれば，(B.4) を満たす．この二次不等式の解の形はやや汚いが，簡便な十分条件とし
て，d − 2 < α ≤ d − 2 + 1/5 が挙げられる．
Appendix C.
定理 3.1 の証明
φ(·; µ, σ 2 ) を d 次元多変量正規分布 Nd (µ, σ 2 Id ) の確率密度関数とする．まず p̂π (z; x)
の分子を θ に関して平方完成すると，w = (nx + z)/(n + 1) に対して
∫
∫
φ(z; θ, 1)φ (x; θ, 1/n) π(θ)dθ = p̂U φ (w; θ, 1/(1 + n)) π(θ)dθ
と書ける．従って，p̂U と p̂π のリスクの差 RKL (θ, p̂U ) − RKL (θ, p̂π ) は
RKL (θ, p̂U ) − RKL (θ, p̂π )
∫
∫
∫
φ (w; θ, 1/(1 + n)) π(θ)dθ
∫
dzdx
= φ (x; θ, 1/n) φ(z; θ, 1) log
φ (x; θ, 1/n) π(θ)dθ
と書ける．さらに w が N (θ, {1/(n + 1)}Id ) に従うことに注意すると，v ∼ N (θ, `Id ) に対
する事前分布 π(θ) の周辺尤度を mπ (v; `) として
RKL (θ, p̂U ) − RKL (θ, p̂π ) = Ew [log mπ (w; 1/{n + 1})] − Ex [log mπ (x; 1/n)]
)
∫ 1/n (
∂
=
− Ev [log mπ (v; `)] d`
∂`
1/(n+1)
(C.1)
168
日本統計学会誌第45巻第1号 2015
と表せる．以下では (C.1) の右辺の被積分関数を書き直す．まず
[
]
∂
∂
∂
Ev [log mπ (v; `)] = Ev
log mπ (v; `) + log mπ (v; `) log φ(v; θ, `)
∂`
∂`
∂`
[
(
)]
kv − θk2
∂
d
= Ew
log mπ (v; `) + log mπ (v; `) − +
∂`
2`
2`2
(C.2)
である．さらに
(
)
(
)
∫
∂
vi − µi
kv − θk2
1
mπ (v; `) =
−
exp
−
π(θ)dθ
∂vi
`
2`
(2π)d/2 `d/2
(
)
∫
∂2
(vi − µi )2
kv − θk2
mπ (v; `)
1
+
exp
−
π(θ)dθ
m
(v;
`)
=
−
π
∂vi2
`
`2
2`
(2π)d/2 `d/2
であり，i = 1, . . . , d をまとめて
d
∇T ∇mπ (v; `) = − mπ (v; `)
`
(
)
∫
1
kv − θk2
kv − θk2
+
exp −
π(θ)dθ
`2
2`
(2π)d/2 `d/2
∂
= 2mπ (v; `) log mπ (v; `)
∂`
(C.3)
を得る．また，スタイン等式 (2.2) を vi ∼ N (θi , `) に注意し，E[kv − θk2 h(v)] に対してく
り返し適用すると，
E[(vi − θi )2 h(v)] = E[(vi − θi ){(vi − θi )h(v)}] = `E[{∂/∂vi }{(vi − µi )h(v)}]
= `E[h(v) + (vi − µi ){∂/∂vi }h(v)] = E[`h(v) + `2 {∂ 2 /∂vi2 }h(v)]
であり，i = 1, . . . , d をまとめて
[
]
E[kv − θk2 h(v)] = E `dh(v) + `2 ∇T ∇h(v)
(C.4)
を得る．よって (C.2) の右辺に (C.3) 及び h(v) = log mπ (v; `) に対する (C.4) を代入す
ると，
[ T
]
∂
1
∇ ∇mπ (v; `)
T
Ev [log mπ (v; `)] = Ev
+ ∇ ∇ log mπ (v; `)
∂`
2
mπ (v; `)
[
]
√
∇T ∇ mπ (v; `)
√
= 2Ev
mπ (v; `)
を得る．
一方 v ∼ Nd (θ, `Id )，尺度について基準化された損失関数 kδ − θk2 /` のもとで，推定量
v と事前分布 π(θ) に対するベイズ推定量 δπ のリスクの差は，(2.15) より
[
]
√
∇T ∇ mπ (v; `)
Ev [kv − θk2 ] − Ev [kδπ − θk2 ]
√
= −4`Ev
`
mπ (v; `)
169
スタイン・パラドクスとベイズ理論
である．従って，
∫
1/n
RKL (θ, p̂U ) − RKL (θ, p̂π ) =
1/(n+1)
(
Ev [kv − θk2 ] − Ev [kδπ − θk2 ]
`
)
d`
2`
を得る．
謝辞
執筆の機会を与えて下さった九州大学の大西先生に感謝する．東京大学の松田孟留さんと
矢野恵佑さんには，有益なコメントを頂いた．ここに記して感謝する．また本研究は JSPS
科研費 25330035 の助成を受けたものである．
参考文献
Aitchison, J. (1975). Goodness of prediction ﬁt, Biometrika, 62(3), 547–554.
Anderson, T. W. (2003). An Introduction to Multivariate Statistical Analysis, 3rd edition, pp. 721, Wiley Series
in Probability and Statistics, Hoboken, NJ, Wiley-Interscience, John Wiley & Sons.
Blyth, C. R. (1951). On minimax statistical decision procedures and their admissibility, Ann. Math. Statist.,
22, 22–42.
Brewster, J. F. and Zidek, J. V. (1974). Improving on equivariant estimators, Ann. Statist., 2, 21–38.
Brown, L. D. (1971). Admissible estimators, recurrent diﬀusions, and insoluble boundary value problems, Ann.
Math. Statist., 42, 855–903.
Brown, L. D. and Hwang, J. T. (1982). A uniﬁed admissibility proof, in Statistical Decision Theory and Related
Topics, III, Vol. 1 (West Lafayette, Ind., 1981 ), pp. 205–230, New York, Academic Press.
Brown, L. D., George, E. I. and Xu, X. (2008). Admissible predictive density estimation, Ann. Statist., 36(3),
1156–1170.
George, E. I., Liang, F. and Xu, X. (2006). Improved minimax predictive densities under Kullback-Leibler loss,
Ann. Statist., 34(1), 78–91.
Ghosh, M. (1994). On some Bayesian solutions of the Neyman-Scott problem, in Statistical Decision Theory
and Related Topics, V (West Lafayette, IN, 1992 ), pp. 267–276, Springer, New York.
James, W. and Stein, C. (1961). Estimation with quadratic loss, in Proc. 4th Berkeley Sympos. Math. Statist.
and Prob., Vol. I , pp. 361–379, Berkeley, Calif., Univ. California Press.
Komaki, F. (2001). A shrinkage predictive distribution for multivariate normal observables, Biometrika, 88(3),
859–864.
久保川達也 (2004). 「スタインのパラドクスと縮小推定の世界」下平英寿，久保川達也，竹内啓，伊藤秀一（編）
『モデル選択予測・検定・推定の交差点』第 3 巻，統計科学のフロンティア，岩波書店，第 3 章．
Lehmann, E. L. and Casella, G. (1998). Theory of Point Estimation, 2nd edition, pp. 589, Springer Texts in
Statistics, New York, Springer-Verlag.
Maruyama, Y. (1998). A uniﬁed and broadened class of admissible minimax estimators of a multivariate normal
mean, J. Multivar. Anal., 64(2), 196–205.
Maruyama, Y. (2003). A robust generalized Bayes estimator improving on the James-Stein estimator for spherically symmetric distributions, Statist. Decisions, 21(1), 69–77.
Maruyama, Y. (2009). An admissibility proof using an adaptive sequence of smoother proper priors approaching
the target improper prior, J. Multivar. Anal., 100(8), 1845–1853.
Maruyama, Y. and Iwasaki, K. (2005). Sensitivity of minimaxity and admissibility in the estimation of a positive
normal mean, Ann. Inst. Statist. Math., 57(1), 145–156.
Maruyama, Y. and Strawderman, W. E. (2005). A new class of generalized Bayes minimax ridge regression
estimators, Ann. Statist., 33(4), 1753–1770.
170
日本統計学会誌第45巻第1号 2015
Maruyama, Y. and Strawderman, W. E. (2013). Improved robust Bayes estimators of the error variance in
linear models, J. Statist. Plann. Inference, 143(6), 1091–1097.
丸山祐造，竹村彰通 (2007). 「統計と特殊関数」『数学セミナー』546, 24–27.
Maruyama, Y. and Takemura, A. (2008). Admissibility and minimaxity of generalized Bayes estimators for
spherically symmetric family, J. Multivar. Anal., 99(1), 50–73.
Stein, C. (1956). Inadmissibility of the usual estimator for the mean of a multivariate normal distribution, in
Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, 1954–1955, vol. I ,
pp. 197–206, Berkeley and Los Angeles, University of California Press.
Stein, C. (1964). Inadmissibility of the usual estimator for the variance of a normal distribution with unknown
mean, Ann. Inst. Statist. Math., 16, 155–160.
Stein, C. (1981). Estimation of the mean of a multivariate normal distribution, Ann. Statist., 9(6), 1135–1151.
Strawderman, W. E. (1971). Proper Bayes minimax estimators of the multivariate normal mean, Ann. Math.
Statist., 42(1), 385–388.

Download Report