Stein’s approach for constructing improved estimators for matrix estimation problems.

今野良彦^{1} 千葉大学理学部数学 · 情報数理学科

0. はじめに多変量正規分布の平均ベクトルの最尤推定量がいわゆる縮小型の推定量で平均２乗誤差

の和の下で一様に改良できること（いわゆる Stein 現象または Stein 効果）が示されて以来、様々な同時推定問題で同様の現象が研究されている . この報告では Stein 現象の中でも matrix estimation problems と呼ばれる推定問題に焦点を絞ってその議論をおこなう . Stein 現象のいろいろな方向への最近の発展に関しては、久保川 (1996)、 Lehmann and Casella (1998)、 Robert (1994), 篠崎 (1991) 等を参照のこと . ここで例として取り上げる推定問題は簡単な推定量を除いて危険関数を直接評価するのが非常に困難である . しかし、推定量の危険関数の直接評価をするのではなく、ある族に含まれる推定量の危険関数にたいする不偏推定量を導出し、個々の推定量の危険関数の不偏推定量を比較することにより、推定量の危険関数を評価する方針が Stein (1973) に提示されてから、正規分布の分散共分散行列の推定問題を含むいくつかの推定問題においてこの方針に沿い、改良型の推定量が提案されている . Stein の方針に沿い、いくつかの推定問題において縮小型の推定量の構成を議論する . 1 節においては、基本的な例である多変量正規分布の平均行列と分散共分散行列の推定問題を定式化し、 2 節では Stein 方針を解説し、 3 節ではこの方針を実行するための技術的な道具をあげ、 4 節では基本的な推定問題をふくむ４つの推定問題において縮小型の推定量の構成を議論する .

1. 推定問題 X : m × p と S : p × p は独立とし、

X ∼ N_{m×p}(μ, T ⊗ Σ),  S ∼ W_{p}(n, Σ)  (1)

とする . ただし、 μ は m × p の未知の平均行列とし、 T は m × m の既知の正定値行列、 Σ は p × p の未知の正定値行列とし、 X の確率密度関数は

\frac{1}{(\sqrt{2π})^{mp}(det T)^{(p/2)}(det Σ)^{(m/2)}} etr [−\frac{1}{2}T^{−1}(x − μ)Σ^{−1}(x − μ)′]

で与えられる . W_{p}(n, Σ) は自由度 n、平均 nΣ の Wishart 分布とする . ここで、 exp tr(·) を etr(·)、行列 A の転置行列を A′ とそれぞれ記した .

未知の母数 μ と Σ の推定問題を個々に考えることにする .

1 email address: konno@math.s.chiba-u.ac.jp この原稿は春季数学会（於：学習院大学）の 1999 年 3 月 28 日のための講演の予稿を訂正したものである .

（訂正日 : 1999 年 3 月 29 日）

例 1 （平均行列の推定）： μ の推定問題を扱うために損失関数として

L_{1}(μ̂, μ; Σ) = tr[Q^{−1}(μ̂ − μ)Σ^{−1}(μ̂ − μ)′]

を採用する . ただし、 Q は既知の m × m の正定値行列で、 μ̂ は (X, S) に基づく μ の推定量である . 例 2 （分散共分散行列の推定）： Σ の推定問題に対して

L_{2}(Σ̂, Σ) = tr(Σ̂Σ^{−1}) − log det(Σ̂Σ^{−1}) − p  (2)

を用いる . ここで、 Σ̂ は Σ の推定量とする .

それぞれの損失関数を (1) に関して期待値を取ったものを危険関数とよび、 R_{1}(μ̂, μ) と R_{2}(Σ̂, Σ) と記すことにする .

注意 1 ：ここでは、推定問題に対してある種の不変性をもつ損失関数のみを考える . 他の損失関数の下での議論として、例 1 では、 Honda (1991)、例 2 では、 Haff (1991) 等がある .

2. 基本的なアプローチ改良型の推定量を導出する基本的な方針は以下の段階を踏むことになる .

1. 考察の対象の統計モデル（尤度関数）を分解し、条件付き独立性等を用いてモデルの分解に対応して推定量や危険関数の分解を行い、より簡単な統計モデルにおける推定問題に帰着させる .

2. モデルの構造を反映した共変性（推移的に作用していないもの！）を用いて、推定量の族を制限する .

3. 共変推定量の族を陽に求める .

4. Normal identity や Wishart identity 等の部分積分の式と calculus on eigenstructure を用いて、求めた共変推定量に対する危険関数の不偏推定量を計算する .

5. 危険関数の不偏推定量から導出された偏微分不等式を利用して、改良型推定量を導く .

注意 2 ：個々の問題でこの方針に沿って改良型推定量を導出するのはモデルの構造に対する深い洞察もしくは直観的な解釈が必要となる . 通常の推定量（たとえば、最尤推定量）は、母数空間に推移的に作用する群のもとでも共変なので、この共変構造に対して不変な損失関数（この報告で採用している損失関数）を採用しているがぎりは定数の危険関数をもつ . しかし、段階 2 において得られる推定量の族は、母数空間に推移的に作用しない群のもとの共変推定量なので、推定量の危険関数を母数に依存するものとなる . また、推定量を制限をするときに、推定量の族を広くとれば、それだけ段階 4 で得られる危険関数の不偏推定量が複雑になり、逆にその扱いが困難になる .

注意 3 ：前章の基本的な推定問題では、段階 1 は不要になる .

注意 4 ：正規分布の平均ベクトルの推定では、危険関数の不偏推定量の評価を経由せずに危険関数をより直接的に評価することは可能であるが、この報告で扱う正規分布の分散共分散行列の推定のように行列母数の同時推定問題では、危険関数の不偏推定量を用いての評価を経由しない直接的な評価は困難である . しかし、例 1 の推定問題において、 Zidek (1978) はある共変推定量の族に対する危険関数を zonal 多項式を用いた表現を与えているが、改良型の推定量を導出するには見通しの良い表現ではなかった .

3. 技術的な補題 \underline{Normal and Wishart identities} . 危険関数の不偏推定量を導出するときに重要な道具が部分積分の式である . この式の意味を理解するために１次元の分布にたいする式をあげる .

補題 1 ：確率変数 Y_{1} は N(ξ_{1}, 1) に従うとし、 g_{1} : ℜ ↦ ℜ は任意の有界閉区間で絶対連続で、 E_{ξ_{1}}[|ġ_{1}(Y_{1})|] < ∞ なる関数とする . このとき、

E_{ξ_{1}}[g_{1}(Y_{1})(Y_{1} − ξ_{1})] = E_{ξ_{1}}[ġ_{1}(Y_{1})]

が成立する . ここで、 ġ_{1} は g_{1} の導関数とした .

補題 2 ：確率変数 V が σ^{2}χ_{n} に従い、 h_{1} : ℜ^{+} ↦ ℜ は任意の有界閉区間で絶対連続で、 E[V^{−1}h_{1}(V)] と E[|ḣ_{1}(V)|] は有限とする . このとき、

E[\frac{h_{1}(V)}{σ^{2}}] = 2E[ḣ_{1}(V)] + (n − 2)E[\frac{h_{1}(V)}{V}]

が成立する . ここで、 χ_{n} は自由度 n のカイ二乗分布とし、 ḣ_{1} は h_{1} の導関数とした .

１次元の正規分布とカイ二乗分布に関する部分積分の式をそれぞれ多次元の分布である多変量変量の正規分布 N(ξ, Σ) と Wishart 分布に関する部分積分の式に拡張したものが Normal identity と Wishart identity である . 前者は Stein (1973) により導入され、後者が Stein (1975) と Haff (1977) に独立に示されている . Normal identity ： Y = (Y_{1}, Y_{2}, …, Y_{p})′ ∼ N_{p}(ξ, Σ) とし、 g : p × 1 = (g_{1}, …, g_{p})′、 g_{i} : ℜ^{p} ↦ ℜ(i = 1, 2, …, p) とする . このとき、

E[Σ^{−1}(Y − ξ)g′(Y)] = E[∇_{Y} g′(Y)]

が成立する . ただし、 ∇_{Y} = (∂/∂Y_{1}, ∂/∂Y_{2}, …, ∂/∂Y_{p})′ である . Wishart identity ： n > p + 1 とする . h : p × p = (h_{ij}) 、 h_{ij} は正定値行列上の実数値関数とし、 S : p × p ∼ W_{p}(n, Σ) とする . このとき、

E[tr{Σ^{−1}h(S)}] = 2E[tr{D_{S}h(S)}] + (n − p − 1)E[tr{S^{−1}h(S)}]

が成立する . ただし、 D_{S} = (d_{ij}) で d_{ij} = (1/2)(1 + δ_{ij})∂/∂S_{ij} である . ここで、 S_{ij} は行列 S の (i, j) 成分とし、

δ_{ij} = \{ \begin{matrix}1 & もし i = j \\ 0 & その他 \end{matrix}

とした . \underline{Calculus on eigenstructure} . S_{1} と S_{2} を p × p の正定値行列とし、 S_{1}S_{2}^{−1} の固有根を f_{1}, f_{2}, …, f_{p}(f_{1} ≥ f_{2} ≥ ⋯ ≥ f_{p} > 0) とかき、 F = diag(f_{1}, f_{2}, …, f_{p}) とおく . A は p × p の正則行列で A′S_{2}A = I_{p} と A′S_{1}A = F を満足するものとする . D_{S_{1}} と D_{S_{2}} を S_{1} と S_{2} の成分をそれぞれ用いて D_{S} と同様に定義したものとする . 更に、実数値関数 φ_{i}(f_{1}, f_{2}, …, f_{p}) は {(f_{1}, f_{2}, …, f_{p}) ∈ ℜ^{p} : f_{1} ≥ f_{2} ≥ ⋯ ≥ f_{p} > 0} 上で定義され、その上で偏微分可能とし、 Φ = diag(φ_{1}, φ_{2}, …, φ_{p}) と φ_{ii} = ∂φ_{i}/∂f_{i}, i = 1, 2, …, p と書くことにする .

補題 3 ：

tr[D_{S_{1}}\{(A′)^{−1}ΦA^{−1}\}] = ∑_{i=1}^{p} \{φ_{ii} + ∑_{j>i} \frac{φ_{i} − φ_{j}}{f_{i} − f_{j}} \}
tr[D_{S_{2}}\{(A′)^{−1}ΦA^{−1}\}] = ∑_{i=1}^{p} \{pφ_{i} − f_{i}φ_{ii} − ∑_{j>i} \frac{f_{i}φ_{i} − f_{j}φ_{j}}{f_{i} − f_{j}} \}
D_{S_{1}}f_{i} = a_{i}a′_{i} i = 1, 2, …, p

ただし、 a_{i} は A の i 列ベクトルとする .

4. リスクの不偏推定量と推定量の構成ここでは、 2 節にある方針に沿い、 1 節にあげた例を含む推定問題について改良型の推

定量を構成する . 例１（続き）： Bilodeau and Kariya (1989) では、 μ の推定量として mp 個の (X, S) 上の実数値関数によって構成される推定量を考えているが、ここでは対象の族を狭めて、

μ̂(X, S) = X + \frac{1}{2}T \{∇_{X}h(F)\}S  (3)

なる形の推定量を考える . ただし、 X = (x_{ij}) に対して、 ∇_{X} : m × p = (∂/∂x_{ij})、 F = diag(f_{1}, f_{2}, …, f_{m∧p}) とし、 f_{1}, f_{2}, …, f_{m∧p}(f_{1} ≥ f_{2} ≥ ⋯ ≥ f_{m∧p} > 0) は X′(T QT)^{−1}XS^{−1} の固有根で、 h は {(f_{1}, f_{2}, …, f_{m∧p}) ∈ ℜ^{m∧p} : f_{1} ≥ f_{2} ≥ ⋯ ≥ f_{m∧p} > 0} 上で定義されている実数値関数で、２回偏微分可能とする . ただし、 m ∧ p = min(m, p) である .

推定量 (3) の危険関数の不偏推定量を表現するために以下を定義する . \underline{記号}：自然数 a, b, c と推定量 (3) に対して、

J(h; a, b, c) = ∑_{i=1}^{c} \{2(b − c + 1)h_{i} + 4f_{i}h_{ii} + 4 ∑_{j>i} \frac{f_{i}h_{i} − f_{j}h_{j}}{f_{i} − f_{j}}

+ (a + c − 3)f_{i}h_{i}^{2} − 4f_{i}^{2}h_{i}h_{ii} − 2 ∑_{j>i} \frac{(f_{i}h_{i})^{2} − (f_{j}h_{j})^{2}}{f_{i} − f_{j}} \}

とおく . ただし、 h_{i} = ∂h/∂f_{i}, h_{ii} = ∂^{2}h/∂f_{i}^{2}(i = 1, 2, …, c) である .

定理 1 ： |p − m| > 1 かつ n > p + 1 とする . このとき、推定量 (3) の危険関数の不偏推定量は

R̂(μ̂, μ) = ptr(QT) + J(h; n + {(m − p) ∧ 0}, m ∨ p, m ∧ p)  (4)

で与えられる . すなわち、

R(μ̂, μ) = E[R̂(μ̂, μ)]

である . ここで、 m ∨ p = max(m, p) とした .

上の定理より、推定量 (3) が X を一様に改良するための条件は J ≤ 0 を満たすような h であればよいことがわかる .

たとえば、 h として、

h = −{(m ∨ p) − (m ∧ p) − 1} ∑_{i=1}^{m∧p} log f_{i}

と取れば、簡単な計算により (4) の右辺の２項目が負になることがわかる . これを書き換えるといわゆる Efron - Morris 型の推定量

μ̂(X, S) = \{ \begin{matrix}X − (m − p − 1)Q^{−1}T^{−1}X\{X′(T QT)^{−1}X\}^{−1}S \\ （m > p + 1 の場合） \\ X − (p − m − 1)T \{XS^{−1}X′\}^{−1}X \\ （p > m + 1 の場合） \end{matrix}

となることがわかる .

注意 5 ：定理 1 におけるパラメータの対応は X′XS^{−1} の固有根の同時確率密度関数のパラメータの対応と同じである .

m > p + 1 の場合  p > m + 1 の場合
p  ↦  m
m  ↦  p
n  ↦  n + m − p.

注意 6 ：定理 1 の結果は分散共分散が既知の場合の多変量正規分布の平均行列の推定問題における Stein (1973) の結果の拡張とみなすことができる . また、 Kuriki (1993) は分散共分散行列が既知の場合の多変量正規分布の p × p の平均行列が歪対称な構造をもつときの平均行列の推定問題において、 Stein (1973) の結果に対応するものを得ている .

注意 7 ： Bilodeau and Kariya (1989) により得られた危険関数の不偏推定量と比較すると推定量の族を制限したことにより、改良型の推定量を構成するにはより見通しの良い式になっている . ただし、彼らが危険行列のもとで X を一様に改良する推定量を提案しているが、そのなかのいくつかの推定量はわれわれの族に含まれないので、この結果は適用できない .

例 2 （続き）： James and Stein (1961) は上三角行列の群のもとで共変な推定量の族（母数空間に推移的に作用するので、その危険関数は定数）の中で最良のものを求めている . Hunt–Stein の定理（Bonder and Milnes (1981) を参照）を利用することで、これはミニマックス推定量であり、ミニマックスリスクが

∑_{i=1}^{p}{log(n + p + 1 − 2i) − E[log χ_{n−i+1}^{2}]}

で与えられることを彼らは示している .

更に、 Stein (1975) では、推定量の族として、直交群に関して共変なものを考え、それが

Σ̂ = Odiag(ψ_{1}, ψ_{2}, …, ψ_{p})O′  (5)

の形で与えられること示している . ここで、

S = OLO′, L = diag(l_{1}, l_{2}, …, l_{p}), l_{1} ≥ l_{2} ≥ ⋯ ≥ l_{p}

を S のスペクトル分解とし、 ψ_{i}(i = 1, 2, …, p) は {(l_{1}, l_{2}, …, l_{p}) ∈ ℜ^{p} : l_{1} ≥ l_{2} ≥ ⋯ ≥ l_{p} > 0} 上で偏微分可能な実数値関数関数とする .

推定量 (5) に対する損失関数は

L(Σ̂, Σ) = tr(Σ̂Σ^{−1}) − ∑_{i=1}^{p} log ψ_{i} + log det(Σ) − p

なので、

ρ(Σ̂, Σ) = E[tr(Σ̂Σ^{−1}) − ∑_{i=1}^{p} log ψ_{i}]

として、この ρ を用いて推定量の比較を行えばよいことがわかる . 更に、 Wishart identity と補題 3 を使えば、 Stein (1975) による ρ の不偏推定量

ρ̂(Σ̂) = ∑_{i=1}^{p} [2ψ_{ii} + (n − p − 1)\frac{ψ_{i}}{l_{i}} + 2 ∑_{i>j} \frac{ψ_{i} − ψ_{j}}{l_{i} − l_{j}} − log ψ_{i}]  (6)

を得る . ただし、 ψ_{ii} = ∂ψ_{i}/∂l_{i}, i = 1, 2, …, p である .

注意 8 ：危険関数の不偏推定量の導出についての証明は竹村 (1991) を参照のこと . また、 Sheena (1995) は、 S の固有根の同時確率密度関数から (6) をより直接的な方法で導出している .

ここで、 (5) において、

ψ_{i} = d_{i}l_{i}, d_{i} = \frac{1}{n + p + 1 − 2i} (i = 1, 2, …, p)  (7)

とおけば、ミニマックス推定量を得る .

しかし、これは ψ_{1} ≥ ψ_{2} ≥ ⋯ ≥ ψ_{p} を常に満足するとはかぎらないので、非許容的になることを Sheena and Takemura (1992) は示し、 (7) で与えられるミニマックス推定量を一様に改良するものを得らている . 一方、 Perron (1992) は上三角行列の群の下で共変なミニマックス推定量を直交行列の群上の一様分布で平均を取り修正をおこなった推定量の近似解を求めることにより、別のタイプのミニマックス推定量を求めている .

もっとも簡単な p = 2 の場合を考えると、 Perron の推定量は (5) において

ψ_{1} = (\frac{l_{1}}{l_{1} + l_{2}} d_{1} + \frac{l_{2}}{l_{1} + l_{2}} d_{2}) l_{1}
ψ_{2} = (\frac{l_{2}}{l_{1} + l_{2}} d_{1} + \frac{l_{1}}{l_{1} + l_{2}} d_{2}) l_{2}

とおけばよいことがわかる . p = 2 の場合において、

ψ_{1} = (\frac{\sqrt{l_{1}}}{\sqrt{l_{1}} + \sqrt{l_{2}}} d_{1} + \frac{\sqrt{l_{2}}}{\sqrt{l_{1}} + \sqrt{l_{2}}} d_{2}) l_{1}
ψ_{2} = (\frac{\sqrt{l_{2}}}{\sqrt{l_{1}} + \sqrt{l_{2}}} d_{1} + \frac{\sqrt{l_{1}}}{\sqrt{l_{1}} + \sqrt{l_{2}}} d_{2}) l_{2}

とすれば、目標の推定量を明示的に与えることを Takemura (1984) が示している . 例 3 （MANOVA – GMANOVA 混合モデルの回帰係数行列の推定）：ここではモデル

Z : n × p = A_{1}β_{1}C_{1} + A_{2}β_{2} + ϵ, ϵ ∼ N_{n×p}(0, I_{n} ⊗ Σ̃)

を考える . ただし、 A_{i} : n × a_{i}(i = 1, 2) と C_{1} は既知の行列でそれぞれのランクを a_{i} と c とする . β_{1} : a_{1} × c_{1} と β_{2} : a_{2} × p は未知の回帰係数行列とする . 更に仮定として、 A : n × (a_{1} + a_{2}) = [A_{1}, A_{2}] のランクは a_{1} + a_{2} とし、 Σ̃ : p × p は未知の正定値行列とする .

(β_{1}, β_{2}) の同時推定問題を考えるために損失関数として、

L((β̂_{1}, β̂_{2}), (β_{1}, β_{2})) = tr[{A_{1}(β̂_{1} − β_{1})C_{1} + A_{2}(β̂_{2} − β_{2})}Σ̃^{−1}{⋯}′]

を採用する .

一方、 (β_{1}, β_{2}) の最尤推定量は

β̂_{1}^{mle} = (A′_{1}P A_{1})^{−1}A′_{1}P ZS̃^{−1}C′_{1}(C_{1}S̃^{−1}C′_{1})^{−1}
β̂_{2}^{mle} = (A′_{2}A_{2})^{−1}A′_{2}[Z − A_{1}β̂_{1}^{mle}C_{1}]

で与えられる . ただし、

S̃ : p × p = Z′[I_{n} − A(A′A)^{−1}A′]Z, P : n × n = I_{n} − A_{2}(A′_{2}A_{2})^{−1}A′_{2}

とした .

推定問題を分解するためにつぎの記号を定義する .

V : a_{1} × p = (A′_{1}P A_{1})^{−1}A′_{1}P Z
U : a_{2} × p = (A′_{2}A_{2})^{−1}A′_{2}Z

としたとき、 V, U, S̃ は互いに独立で

V ∼ N_{a_{1}×p}(β_{1}C_{1}, (A′_{1}P A_{1})^{−1} ⊗ Σ̃)  (8)
U ∼ N_{a_{2}×p}(η, (A′_{2}A_{2})^{−1} ⊗ Σ̃)  (9)
ただし、 η = β_{2} + (A′_{2}A_{2})^{−1}A′_{2}A_{1}β_{1}C_{1}
S̃ ∼ W_{p}(n − (a_{1} + a_{2}), Σ̃)

に従う . 更に、 G : c_{1} × p = (C_{1}C′_{1})^{−1}C_{1} とおき、 G_{o} : (p − c_{1}) × p は GG′_{o} = 0 かつ G_{o}G′_{o} = I_{p−c_{1}} を満足するものとし、

V_{1} : a_{1} × c_{1} = V G′, V_{2} : a_{2} × (p − c_{1}) = V G′_{o}

とおき、

S = [ \begin{matrix}S_{11} & S_{12} \\ S_{21} & S_{22} \end{matrix} ] = [ \begin{matrix}GS̃G′ & GS̃G′_{o} \\ G_{o}S̃G′ & G_{o}S̃G′_{o} \end{matrix} ]

とし、

S_{11·2} = S_{11} − S_{12}S_{22}^{−1}S_{21}

と定める . また、 Σ と Σ_{11·2} も G と G_{o} を用いて Σ̃ から同様に定める .

(V_{1}, V_{2}, S_{11·2}, S_{22}^{−1}S_{21}, S_{22}) と (U, S) はそれぞれ (β_{1}, Σ) と (β_{2}, Σ) の十分統計量なので、推定量の族として

β̂_{1}(B_{1}) = β̂_{1}^{mle} + B_{1}(V_{1}, V_{2}, S_{11·2}, S_{22}^{−1}S_{21}, S_{22})
β̂_{2}(B_{2}) = β̂_{2}^{mle} + B_{2}(U, S)

を考える . ただし、 B_{1} : a_{1} × c_{1} と B_{2} : a_{2} : ×p とする . すると β̂_{1}(B_{1}) と β̂_{2}(B_{2}) の危険関数は

R((β̂_{1}(B_{1}), β̂_{2}(B_{2})), (β_{1}, β_{2})) = R_{1} + R_{2}

となる . ここで、

R_{1} = E[tr{(A′_{1}P A_{1})(β̂_{1}^{mle} + B_{1} − β_{1})Σ_{11·2}^{−1}(⋯)′}]
R_{2} = E[tr{(A′_{2}A_{2})(U + B_{2} − η)Σ^{−1}(U + B_{2} − η)′}]

となる . これより、 (β_{1}, β_{2}) の同時推定問題はふたつの独立した推定問題に分解ができる .

[V_{1}, V_{2}] ∼ N_{a_{1}×p}([β_{1}, 0], (A′_{1}P A_{1})^{−1} ⊗ Σ)  (10)

と (9) に注意すれば、ひとつは GMANOVA モデル (10) における回帰係数行列 β_{1} の推定を危険関数 R_{1} のもとで考えた問題となり、他方は MANOVA モデル (9) における回帰係数行列 η の推定を危険関数 R_{2} のもとの問題となる . 後者については例 1 の結果を直接適応可能であり、前者については GMANOVA モデルを条件付き独立性で分解し、それに対応する推定量の族を考えることで β_{1} の推定問題をさらに分解し、 β_{1} + V_{2}Σ_{22}^{−1}Σ_{21} と Σ_{22}^{−1}Σ_{21} に関する (V_{2}, S_{22}) を条件付きした推定問題として別個に扱うことができることを Kariya, et al (1996) は示している . 例 4 （LCI 構造をもつ分散共分散行列の推定）：集合 I = {1, 2, …, p} に対して、 K を I の部分集合からなる分配環とし、（集合の）和と積について閉じて空集合を含むものとする . X ∼ N_{p}(0, Σ) とし、 K ∈ K に対して、 X_{K} を K の元に対応する座標の成分をもつ部分ベクトルとする . このとき、 Lattice conditional independence ( LCI ) モデル N(K) を、

すべての L, M ∈ K に対して、 X_{L∩M} を与えたときに
X_{L} と X_{M} は条件付き独立  (11)

であるような多変量正規分布 N(0, Σ) の集合とする . また、 P_{I}(K) をすべての L, M ∈ K に対して、 (11) を満足するような多変量正規分布の分散共分散行列の集まりとする .

LCI モデル N(K) からのランダム標本 X_{1}, X_{2}, …, X_{n}(n ≥ p) に基づいて分散共分散行列 Σ の推定問題を損失関数 (2) のもとで考える . \underline{記号}： K ∈ K に対して、 ⟨K⟩ = ∪(K̃ ∈ K|K̃ ⊂ K かつ K̃ ≠ K) 、 [K] = K\⟨K⟩ とする . したがって、 K = ⟨K⟩ ∪ [K] かつ ⟨K⟩ ∩ [K] は空集合である . 更に、

J (K) = {K ∈ K|[K] は空でない}

とおく . また、 M_{I} と P_{I} を p × p の正方行列と正定値行列の集合とする . また、 A = (a_{i_{1}i_{2}}|(i_{1}, i_{2}) ∈ I × I) ∈ M_{I} と K, L ∈ I に対して、

A_{K×L} = (a_{i_{1}i_{2}}|(i_{1}, i_{2}) ∈ K × L)

とし、その集合を M_{K×L} と書く . また、 A_{K×L} と M_{K×K} を A_{K} と M_{K} とそれぞれ書くことにし、 P ∈ P_{I} に対して同様に P_{K} と P_{K} を定義する .

A ∈ M_{I}、 K ∈ K とし、 K = [K] ∪ ⟨K⟩ に対応して A_{K} の分割を

( \begin{matrix}A_{⟨K⟩} & A_{⟨K]} \\ A_{[K⟩} & A_{[K]} \end{matrix} )

とし、 A_{⟨K⟩} ∈ M_{⟨K⟩} 、 A_{[K]} ∈ M_{[K]} 、 A_{⟨K]} ∈ M_{⟨K⟩×[K]} 、 A_{[K⟩} ∈ M_{[K]×⟨K⟩} とした . また、 Σ ∈ P_{I} と K ∈ J (K) に対して、 Σ_{[K]•} = Σ_{[K]} − Σ_{[K⟩}Σ_{⟨K⟩}^{−1}Σ_{⟨K]} とする .

行列 A ∈ M_{I} が K - preserving であるとは、すべての K ∈ K と x ∈ ℜ^{p} に対して、 (Ax)_{K} = A_{K}x_{K} が成立することである . これは

L, M ∈ J (K) : M ⊈ L ならば A_{[L]×[M]} = 0

が成立することと同値である . GL_{I}(K) と GT_{I}(K) をそれぞれ K - preserving で正則な p × p 行列と正の値の対角成分をもつ K - preserving な下三角行列の群とする .

つぎの性質は推測問題の議論の上で基本的な役割を果たす .

(i) P_{I}(K) ↦ X(M_{[K]×⟨K⟩} × P_{[K]}|K ∈ J (K)) は bijective map である .

(ii) Σ ∈ P_{I}(K) であるための必要十分条件は、 x ∈ ℜ^{p} に対して、

tr(Σ^{−1}xx′) = ∑(tr{Σ_{[K]•}^{−1}(x_{[k]} − Σ_{[K⟩}Σ_{⟨K⟩}^{−1}x_{⟨K⟩})(⋯)′}|K ∈ J (K))

を満足することである .

(iii) Σ ∈ P_{I}(K) と L ∈ K に対して、

det(Σ_{L}) = Π(det Σ_{[K]•}|K ∈ J (K), K ⊆ L)

が成立する .

(iv) GL_{I}(K) による P_{I}(K) への作用 (A, Σ) ↦ AΣA′ は推移的である . ただし、 A ∈ GL_{I}(K), Σ ∈ P_{I}(K) である .

S = ∑_{i=1}^{n} X_{i}X′_{i} とおく . Σ ∈ P_{I}(K) から LCI モデルの性質 (i) を使い、 Σ の最尤推定量は

nΣ̂_{[K]•} = S_{[K]•}, Σ̂_{[K⟩}Σ_{⟨K⟩}^{−1} = S_{[K⟩}S_{⟨K⟩}^{−1}, K ∈ J (K)  (12)

によって一意的に構成できることを Andersson and Perlman (1993) が示している . (12) によって得られる Σ の最尤推定量を Σ̂_{K}^{mle} ∈ P_{I}(K) と記すことにする . 例 2 の場合と同様、 GT_{I}(K)–共変な推定量の中で最良な推定量がミニマックス推定量になることがわかる . そこで、 nΣ̂_{K}^{mle} = T T′ と分解する . ただし、 T ∈ GT_{I}(K) である . 更に、 D = diag(D_{[K]}|K ∈ J (K)) で D_{[K]} は |[K]| × |[K]| の対角行列でその対角成分は

d_{[K],i}^{−1} = n + |[K]| − |⟨K⟩| + ∑(|[L]||L ∈ M(K)) − 2i + 1,

である . ただし、 i = 1, …, |[K]|, K ∈ J (K) 、 M(K) = {L ∈ J (K)|K ⊆ ⟨L⟩} とし、 K ∈ K に対して、 |K| を K の中に含まれる要素の個数とする .

定理 2 ： Σ̂_{K}^{m} = T DT′ はミニマックスである . そのミニマックスリスクは

∑_{K} { ∑_{i=1}^{|[K]|} (− log d_{[K],i} − E log χ_{n+|[K]|−|⟨K⟩|−2i+1}^{2}) |K ∈ J (K)}

で与えられる .

更に、 max{|[K]||K ∈ J (K)} ≥ 2 の場合には Σ̂_{K}^{m} は一様に改良できる . S_{[K]•} をスペクトル分解して、

S_{[K]•} : |[K]| × |[K]| = O_{[K]}diag(l_{[K],1}, l_{[K],2,}, …, l_{[K],|[K]|})O′_{[K]}

と書く . ただし、 l_{[K],1} ≥ l_{[K],2,} ≥ ⋯ ≥ l_{[K],|[K]|} > 0 である .

ここで、以下の方程式から P_{I}(K) に含まれることから LCI モデルの性質 (i) を用いて構成される推定量 Σ̂_{K} を考える .

nΣ̂_{[K]•} = O_{[K]}diag(ψ_{[K],1}, ψ_{[K],2}, …, ψ_{[K],|[K]|})O′_{[K]}
Σ̂_{[K⟩}Σ_{⟨K⟩}^{−1} = S_{[K⟩}S_{⟨K⟩}^{−1}, K ∈ J (K).  (13)

ただし、 ψ_{[K],i}(i = 1, 2, …, |[K]|, K ∈ J (K)) は

{(l_{[K],1}, l_{[K],2,}, …, l_{[K],|[K]|}) ∈ ℜ^{|[K]|}|l_{[K],1} ≥ l_{[K],2,} ≥ ⋯ ≥ l_{[K],|[K]|} > 0}

上で定義され、その上で偏微分可能な実数値関数とする .

LCI モデルの性質を利用して、 (13) で与えれる推定量 Σ̂ の危険関数を分解することにより、例 2 の結果が適用できる .

定理 3 ： Σ ∈ P_{I}(K) とする .

ρ̂(Σ̂_{K}, Σ) = ∑_{K}\{ ∑_{i=1}^{|[K]|}(2ψ_{[K],ii} + (e_{[K]} − |[K]| − 1)\frac{ψ_{[K],i}}{l_{[K],i}}

+2 ∑_{j>i} \frac{ψ_{[K],j} − ψ_{[K],i}}{l_{[K],j} − l_{[K],i}} − log ψ_{[K],i})|K ∈ J (K)\}.

ただし、 ψ_{[K],ii} = ∂ψ_{[K],i}/∂l_{[K],i}, i = 1, 2, …, |[K]| 、 e_{[K]} = n − |⟨K⟩| + {∑ |[L]||L ∈ M(K)}, K ∈ J (K) である .

定理 4 ： max{|[K]||K ∈ J (K)} ≥ 2 のとき、

ψ_{[K],i} = d_{[K],i}l_{[K],i}, i = 1, 2, …, |[K]|, K ∈ J (K)

と (13) においておけば、 Σ̂_{K} はミニマックスで Σ̂_{K}^{m} を一様に改良する .

I = {1, 2, 3} の場合についてふたつの LCI モデルにおけるミニマックス推定量を定理 2 を使い求めてみる . \underline{その１}：K_{1} = \{ϕ, 1, 12, 13, I\} とする . したがって、 J (K_{1}) = \{1, 12, 13\} となる . X_{1}, X_{2}, …, X_{n} を N(K_{1}) からのランダム標本とし、 S = ∑_{i=1}^{n} X_{i}X′_{i} = (s_{ij}) とおく . LCI モデルの性質 (i) から Σ の最尤推定量を求めると

nΣ̂‘_{K_{1}}^{mle} = [ \begin{matrix}s_{11} & s_{12} & s_{13} \\ s_{21} & s_{22} & s_{21}s_{11}^{−1}s_{13} \\ s_{31} & s_{31}s_{11}^{−1}s_{12} & s_{33} \end{matrix} ]

となる . これを分解して、 nΣ̂_{K_{1}}^{mle} = T_{1}T′_{1} と書く . ただし、

T_{1} = [ \begin{matrix}s_{11}^{1/2} & 0 & 0 \\ s_{21}s_{11}^{−1/2} & t_{22} & 0 \\ s_{31}s_{11}^{−1/2} & 0 & t_{33} \end{matrix} ] ∈ GT_{I}(K_{1})

である . ここで、 t_{ii}^{2} = s_{ii} − s_{i1}s_{11}^{−1}s_{1i} とおき、 t_{ii} > 0(i = 2, 3) である . ミニマックス推定量は Σ̂_{K_{1}}^{m} = T_{1}D_{1}T′_{1} で与えられる . ただし、

D_{1} = diag (\frac{1}{n + 2}, \frac{1}{n − 1}, \frac{1}{n − 1}) .

である . 更に、これを成分で書けば、

Σ̂_{K_{1}}^{m} = [ \begin{matrix}\frac{s_{11}}{n+2} & \frac{s_{12}}{n+2} & \frac{s_{13}}{n+2} \\ \frac{s_{21}}{n+2} & \frac{1}{n−1}(s_{22} − \frac{3}{n+2}s_{21}s_{11}^{−1}s_{12}) & \frac{s_{21}s_{11}^{−1}s_{13}}{n+2} \\ \frac{s_{31}}{n+2} & \frac{s_{31}s_{11}^{−1}s_{12}}{n+2} & \frac{1}{n−1}(s_{33} − \frac{3}{n+2}s_{31}s_{11}^{−1}s_{13}) \end{matrix} ]

となることがわかる . \underline{その２}： K_{2} = \{ϕ, 1, 2, 12, I\} とする . したがって、 J (K_{2}) = \{1, 2, I\} となる . 同様に X_{1}, X_{2}, …, X_{n} を N(K_{2}) からのランダム標本とする . Σ(∈ P_{I}(K_{2})) の最尤推定量は nΣ̂_{K_{2}}^{mle} = T_{2}T′_{2} で与えられる . ここで、

T_{2} = [ \begin{matrix}s_{11}^{1/2} & 0 & 0 \\ 0 & s_{22}^{1/2} & 0 \\ t_{31} & t_{32} & t_{33} \end{matrix} ] ∈ GT_{I}(K_{2})

とし、

(t_{31}, t_{32}) = (\frac{s_{31}}{\sqrt{s_{11}}} , \frac{s_{32}}{\sqrt{s_{22}}} ), t_{33}^{2} = s_{33} − (\frac{s_{31}^{2}}{s_{11}} + \frac{s_{32}^{2}}{s_{22}} )

とおいた . ミニマックス推定量は Σ̂_{K_{2}}^{m} = T_{2}D_{2}T′_{2} で与えられる . ただし、

D_{2} = diag (\frac{1}{n + 1}, \frac{1}{n + 1}, \frac{1}{n − 2})

である .

5. むすびにかえて 4 節で取り上げた例以外にも matrix estimation problems に関する研究は多々ある . その

中でも特に Stein の着想に沿って改良型の推定量を導出した研究として以下のものがある . 多変量正規分布の２標本問題においる分散共分散行列の同時推定問題（Loh (1988, 1991a, 1991c) ）、同じ問題を別の方向からの定式化したもの（Muirhead and Verathaworn (1985), Bilodeau and Srivastava (1989)）、多変量正規分布の２標本問題における共通平均の推定問題（Loh (1988, 1991b)）、付加情報がある場合の多変量正規分布の共分散行列の推定問題（Konno (1995), Krishinamoorthy (1992)）、 4 節の例 1 と 2 の結果を多変量 t 分布や楕円型分布のもとに拡張した研究（Joarder and Ali (1997), Kubokawa (1997)）、非心 Wishart 分布と非心多変量 F 分布の非心母数の推定問題（Leung (1994a, 1994b) ）等がある .

最近、 Anderson (1984) に代表される多変量解析にあるモデルを広く拡張した多変量統計モデルを提案とその上で推測理論の展開が Andersson and Perlman (1993) の他に Andersson and Madsen (1998)、 Massam and Neher (1998) がある . これらの論文では invariance がモデルの構成に重要な役割を果たしているので、この新しい多変量統計モデルのもとでも Stein の着想はうまく機能することが期待される .

ミニマックス推定という立場から matrix estimation problems の取り組む研究は多いが、 Stein の方針という手段以外に有効に機能する道具がこの問題ではないことから別の方向からの研究はまだまだ数が少ない . 多変量正規分布の分散共分散行列の推定問題において Bayes 推定の観点からは Haff (1991) による the variational form of Bayes estimators という結果や Leonard and Hsu (1992) や Yang and Berger (1994) の論文があるが、許容性という視点からはほとんど研究がされていない .

参考文献

Anderson, T.W. (1984): Introduction to Multivariate Statistical Analysis, 2nd edtion. New York Wiley.

Andersson, S.A. and Madsen, J. (1998): Symmetry and lattice conditional independence in a multivariate normal distribution, Ann. Statist. 26, 525–572.

Andersson, S.A. and Perlman, M.D. (1993): Lattice models for conditional independence in a multivariate normal distribution, Ann. Statist. 21, 1318–1358.

Bilodeau, M., and Kariya, T. (1989): Minimax estimators in the normal MANOVA model, J. Multivariate Anal. 28, 260–270.

Bilodeau, M. and Srivastava, M.S. (1992): Estimation of the eigenvalues of Σ_{1}Σ_{2}^{−1}, J. Multivariate Anal. 41, 1–13.

Bonder, J. and Milnes, P. (1981): Amenability: A survey for statistical applications of Hunt-Stein and related conditions on groups, Z. Wahrsch. verw. Gebiete 57, 103–128.

Dey, D.K. and Srinivasan, C. (1985): Estimation in a covariance matrix under Stein’s loss, Ann. Statist. 13, 1581–1591.

Eaton, M.L. ( 1970): Some problems in covariance estimation. Technical Report No. 49, Department of Statistics, Stanford University.

Efron, B. and Morris, C. (1976): Families of minimax estimators of the mean of a multivariate normal distribution, Ann. Statist. 4, 11–21.

Gleser, L. (1986): Minimax estimators of a normal mean vector of a normal distribution with unknown covariance matrix, Ann. Statist. 14, 1625–1633.

Haff, L. R. (1976): Estimation of the inverse covariance matrix: Random mixtures of the inverse Wishart matrix and identity, Ann. Statist. 6, 1264–1276.

Haff, L. R. (1979): Identity for the Wishart distribution with applications, J. Multivariate Anal. 9, 531–544.

Haff, L. R. (1991): The variational form of certain Bayes estimators, Ann. Statist. 19, 1163–1190.

Honda, T. (1991): Minimax estimators in the MANOVA model for arbitrary quadratic loss and unknown covariance matrix, J. Multivariate Anal. 36, 113–120.

James, W. and Stein, C. (1961): Estimation with quadratic loss, in Proc. Fourth Berkeley Symp. Math. Statist. Probab 1, Univ. California Press, 361–380.

Joarder, A.H. and Ali, M.M. (1997): Estimation of the scale matrix of a multivariate T-model under entropy loss. Metrika 46, 21–32.

Kariya, T., Konno, Y. and Strawderman, W.E. (1996): Double Shrinkage estimation in the GMANOVA model, J. Multivariate Anal. 56, 245–258.

Kariya, T., Konno, Y. and Strawderman, W.E. (1999): Construction of shrinkage estimators for the regression coefficient matrix in the GMANOVA models, to appear in Communication in Statistics, Theory and Method.

Konno, Y. (1991): On estimation of a matrix of normal means with unknown covariance

matrix, J. Multivariate Anal., 36, 44–55. Konno, Y. (1995): Estimation of a normal covariance matrix with incomplete data under

Stein’s loss, J. Multivariate Anal. 52, 308–324. Konno, Y. (1998): Inaddmissibility of the maximum likelihood estimator of normal covari-

ance matrices with the lattice conditional independence, Technical Report No. 340, Department of Statistics, University of Washington.

Konno, Y., Kubokawa, T. and Saleh, A.K.Md.E. (1997): Shrinkage estimators in a mixed

MANOVA and GMANOVA mdoel, Statistics & Decisions 15, 37–49. Krishnamoorthy, K. (1992): Estimation of normal covariance matrix and precision matrices

with incomplete data, Comm. Statist. Theory Methods 20, 757–770. 久保川達也 (1996): 縮小推定の理論と応用 (1), 経済学論集 61, no. 4, 2–31. Kubokawa, T., Saleh, A.K.Md.E., and Morita, K. (1992): Improving on MLE of coefficient

matrix in a growth curve model, J. Statist. Plan. and Infer. 31, 169–177. Kubokawa, T. and Srivatava, M.S. (1997): Robust improvements in estimation of mean

and covariance matrices in ellipitically contoured distribution, 「推定論とその応用の研究」シンポジウム予稿集, 1–10.

Kuriki, S. (1993): Orthogonally invariant estimation of the skew-symmetric normal mean

matrix. Ann. Inst. Statist. Math. 45, 731–739. Lehmann, E. L. and Casella, G. (1998): Theory of Point Estimation, 2nd edtion. Springer. Leonard, T. and Hsu, J.S.J. (1992): Bayesian inference for a covariance matrix.

it Ann. Statist. 20, 1669–1690. Leung, P.-L. (1994a): An identity for the noncentral Wishart distribution with application.

J. Multivariate Anal. 48, 107–114. Leung, P.-L. (1994b): An identity for the noncentral multivariate F distribution with

application. Unpublished manuscript. Loh, W.-L. (1988): Estimating covariance matrices. Ph.D. thesis. Stanford University. Loh, W.-L. (1991a): Estimating covariance matrices, Ann. Statist. 19, 293–296. Loh, W.-L. (1991b): Estimating the common mean of two multivariate normal distribu-

tions, Ann. Statist. 19, 297–323. Loh, W.-L. (1991c): Estimating covariance matrices II, J. Multivariate Anal. 36, 163–174. Massam, H. and Neher, E. (1998): Estimation and testing for lattice conditional indepen-

dence models on Euclidean Jordan algebras, Ann. Statist. 26, 1051–1082. Muirhead, R.J. (1982): Aspects of multivariate statistical analysis, Wiley. Muirhead, R.J. and Verathaworn, T. (1985): Estimating the latent roots of Σ_{1}Σ_{2}^{−1}, In

Multivariate Analysis VI, ( ed. Krishnaiah, P.R.), 431–438. Perron, F. (1992): Minimax estimators of a covariance matrix, J. Multivariate Anal. 43,

6–28. Robert, C.P. (1994): The Bayesian Choice. Springer-Verlag. Sheena, Y. (1995): Unbiased estimator of risk for an orthogonally invariant estimator of a

covariance matrix. J. Japan Statist. Soc. 25, 35–48. Sheena, Y and Takemura, A. (1992): Inadmissibility of non-order-preserving orthogonally

invariant estimators of the covariance matrix in the case of Stein’s loss, J. Multivariate Anal. 41, 117–131.

篠崎信雄 (1991): Stein タイプの縮小推定量とその応用, 応用統計学 20, 59–76. Stein, C. (1956): Inadmissibility of the usual estimator for the mean of a multivariate

distribution. Proc. Third Berkeley Symp. Math. Statist. Prob. 1, University of California Press, 197–206.

Stein, C. (1973): Estimation of mean of a multivariate normal distribution, in Proceedings

of Prague Symposium on Asymptotic Statistics, pp. 345–381. Stein, C. (1975): Estimation of a covariance matrix. Rietz Lecture, 39th Annual Meeting

IMS. Atlanta, Georgia. Takemura, A. (1984): An orthogonally invariant minimax estimator of the covariance

matrix of a multivariate normal population. Tsukuba J. Math. 8, 367–376. 竹村彰通 (1991): 多変量推測統計の基礎, 共立出版 . Tan, M. (1991): Improved estimators for the GMANOVA problem with application to

Monte Carlo Simulation, J. Multivariate Anal. 38, 262–274. Yang, R. and Berger, J.O. (1994): Estimation of a covariance matrix using the refernce

prior. Ann. Statist. 22, 1195–1211. Zheng, Z. (1986): On estimation of matrix of normal mean, J. Multivariate Anal. 18,

70–82.