主成分分析における固有値・固有ベクトルの信頼区間〜漸近理論およびリサンプリング法に基づく信頼区間の性能評価〜
M2007MM016 松下剛士
指導教員：田中豊

1

はじめに

主成分分析において固有値や固有ベクトルは、点推定だけが利用される事が多い。しかし、主成分分析における信頼区間・信頼領域と求める事が出来れば、主成分分析に使われている固有値固有ベクトルの信頼度を知ることができ、主成分をどこまで取り上げるか、また、主成分の解析をどこまで細かく行うかの参考にすることができる。
主成分分析における固有値・固有ベクトルの信頼区間については、多変量正規分布を仮定したときの漸近理論に基づく方法のほか、ブートストラップなどのリサンプリング法に基づく方法が知られている。本研究はそれらの方法による信頼区間の性能をシミュレーションを用いて評価することを目的とする。
実験に用いるデータは、[1] の P84 表 3 のデータから取られる平均と分散共分散行列を用いる。そこで得られた平均と分散共分散行列を母集団の平均と分散共分散行列として、mvrnorm でサンプルデータを作成し、漸近信頼区間やブートストラップ信頼区間を求める。その信頼区間を用いて被覆率を出して、性能評価を行う。
漸近信頼区間には [5] を、ブートストラップ信頼区間には [2][3][6]、t 法のジャックナイフ分散推定量は [4] を参考にした。

2lklk≤ λk ≤1 + τ zα
1 − τ zαと表わされる。この信頼区間を引き出す時に、n が十分に大きいと仮定し、τ zα ≤ 1 となりたつものとする。
別の方法で近似信頼区間を出すには、log lk の分布を用いる方法がある。式 (1) に対してデルタ法を適用すると、

log lk ∼ N (log λk ,2)n−1(3)となる。これによって、分散が未知のパラメータ λk に依存する性質を取り除くことになる。信頼係数 1 − α を持つ log λk のための近似信頼限界は、log lk ± τ zα となり、λk への変換によって、

lk e−τ zα ≤ λk ≤ lk eτ zα(4)という近似信頼区間を得る事ができる。
lk は、漸近的に独立である。よって、λk のいくつかに対する同時信頼領域は、全体的に一定の信頼水準を得るために個々の信頼係数を選び、式 (2) や式 (4) の様な信頼区間を、単純に組み合わせる事によって得る事ができる。

2.2

漸近理論による固有ベクトルの信頼区間

αkj のそれぞれの近似信頼区間は、漸近分散共分散がp∑
λl αlj αlj ′λkn−1(λl − λk )2
l=1,l̸=kcov(akh , ak′ j ′ ) =λk λk′ αkj αk′ j ′−(n − 1)(λk − λk′ )2

漸近信頼区間(2)′k=k
この節における l は多変量正規分布に従う母集団から抽出されたサンプルの分散共分散行列 S の固有値であり、k ̸= k ′a は S の固有ベクトルである。また、λk と αk はそれぞれ、λk = E(lk ) と αk = E(ak ) であり、zα は、標準正規となることを利用して akj の周辺分布から作られる。区分布の上側 100α% 点である。
間は λk の時と同様の方法で作られるが、個々の akj は互
2.1 漸近理論による固有値の信頼区間いに独立でないため、いくつかの αkj のために同時信頼固有値の漸近分散共分散は、領域を見つける必要があり、式はより難くなる。ak の周2辺分布は、
 2λkk = k′cov(lk , lk′ ) = n − 1ak 〜N (αk , T k )0k ̸= k ′と近似することができる。ここで、であり、lk に対する周辺分布は、

2λ2k)lk ∼ N (λk ,n−1(1)と近似することができる。これより、

lk − λ k〜N (0, 1)
λk [2/(n − 1)]1/2を導く事が出来る。τ 2 = 2/(n − 1) とし、ここから λk に対する信頼係数 (1 − 2α) の信頼区間を導くと、

Tk =p∑′
λkλlαl αln−1(λl − λk )2l=1,l̸=kである。行列 T k は、固有ベクトル ak に対応する 1 つのゼロ固有値を持つので、ランクは (p − 1) である。最終的には一般逆行列を用いて、近似的に′(n − 1)αk (lk S −1 + lk−1 S − 2I p )αk ≤ χ2(p−1);α(5)となる。ただし、今回は母集団の分散共分散がわかる
ため、aki ∼ N (αki , T k(i,i) ) と式 (2) よりaki −√√Tk(i,i) zα ≤ αki ≤ aki + Tk(i,i) zα(6)を用いる。Tk(i,i) は ak からなる分散共分散行列 T k の i番目の対角要素である。

3ブートストラップ信頼区間

この節では信頼区間を構成するためのブートストラップ法について述べる。ブートストラップ法を用いた信頼区間のもっとも簡単な方法は、パーセンタイル法である。
次にブートストラップ t 法があるが、この方法では推定量の分散の推定が必要になる。分散があまり信頼できない場合には、この方法を用いるのに注意を要する。最後
図 1 ブートストラップ標本の作成の方法は、BCa 法と呼ばれるもので、パーセンタイル法を改良した方法で、推定量の偏りとその分布の歪みを同求める必要がないので簡単な方法だが、信頼区間の被覆時に補正する。
以下のブートストラップ信頼区間はノンパラメトリッ 率を見ると満足できるものだというわけではない。一般
ク・ブートストラップ法のもとで得られる信頼区間である。 的には本節で述べるブートストラップ t 法のほうが精度がいいと言われている。この節では θ̂ の標準誤差に対す√3.1 信頼区間についてる推定値 σ̂/ n が得られるとする。式 (7) より、スチュー√
信頼区間の正確度について参考文献 [3] より、パーセン デント化された量 T = n(θ̂ − θ)/σ̂ の分布に対しての正規近似は、1 次の正確度しか持っていないが、ブートストタイル法の正確度は、ラップ標本に基づく T ∗ の分布は、T の分布に対し 2 次の√P r{ n(θ̂ − θ)/σ̂ ≤ t} = Φ(t) + O(n−1/2 )(7)正確度を持っている。よってで表わされ、t 法や BCa 法の正確度は、

P r{θ ∈ IL } − (1 − α) = O(n−1 )(θ̂ − n−1/2 σ̂w1−α , θ̂ − n−1/2 σ̂wα )(8)と表わされている。
式 (7) では、1 次の正確度を持っているといい、式 (8)では、2 次の正確度を持っているという。
パーセンタイル法

3.2

パーセントタイル法による信頼区間は、ブートストラップサンプルを小さい順に並び変えて、求めたい点にある値を信頼限界とするものである。実際には、任意の α に対し、ブートストラップ反復回数を B とすると、ブートストラップサンプルの Bα 番目の点を求めればいい。

3.2.1

アルゴリズム(i). 図 1 の様に、もとの標本 y1 , . . . , yn から重複を許して無作為に大きさ n の標本 y1∗ , . . . , yn∗ を抽出をする事で作成する。
(ii). ブートストラップ推定量 θ̂∗ = θ(Y1∗ , . . . , Yn∗ ) の値の∗計算を B 回繰り返し、θ̂1∗ , . . . , θ̂Bを求める。
∗∗とする。Bα は自然数とし、θ̂(Bα)
は(iii). θ̂α = θ̂(Bα)∗θ̂1∗ , . . . , θ̂Bに対する Bα 番目の順序統計量とする。
このとき (θ̂α , θ̂1−α ) は、θ に対する名目上の被覆率1 − 2α をもつ両側信頼区間とする。

3.3

ブートストラップ t 法

上で述べたパーセンタイル法は、ブートストラップ信頼区間を構成するために用いる推定量 θ̂ の分散推定値をは 2 次の正確度を持つ信頼度 1 − 2α の信頼区間になる。

3.3.1

ジャックナイフ分散推定

上の t 法で固有値と固有ベクトルの分散を求める時、通常の方法では分散を求める事が出来ない。そこで、ここでは [4] にあるジャックナイフ分散推定量を用いる。
(i)θ̃n,i = nθ̂n − (n − 1)θ̂n−1 ,i = 1, . . . , n(i)をジャックナイフ疑似量と言う。θ̂n−1 は、元のデータから i 番目のデータを抜いて計算した θ̂ である。疑似量θ̃n,i , i = 1, . . . , n は互いに独立で同じ分布に従う確率変数√のように扱える、θ̃n,i の分散は nθ̂n の分散と近似的に等√しい、という事が正しいとして提案された分散 V ( nθ̂n )の推定量
nn1 ∑1∑(θ̃n,i −θ̃n,j )2n−1ni=1j=1√
より、V (θ̂n ) = V ( nθ̂n )/n の推定量はnn∑
11∑(θ̃n,i −θ̃n,j )2n(n − 1) i=1n j=1(9)nnn − 1 ∑ (i)1 ∑ (j) 2=(θ̂θ̂−)
n i=1 n−1 n j=1 n−1(10)

VJACK =

3.3.2

アルゴリズム3.4.1(i). 図 1 の様に、もとの標本 y1 , . . . , yn から重複を許して無作為に大きさ n の標本 y1∗ , . . . , yn∗ を抽出をする事で作成する。

アルゴリズム(i). 図 1 の様に、もとの標本 y1 , . . . , yn から重複を許して無作為に大きさ n の標本 y1∗ , . . . , yn∗ を抽出をする事で作成する。

(ii). (i) で得られたブートストラップ標本を使い、θ̂ と σ̂ (ii). ブートストラップ推定量 θ̂∗ = θ(Y1∗ , . . . , Yn∗ ) の値の∗を求める。に対応するブートストラップ推定値 θ̂∗ 、σ̂ ∗ を計算
計算を B 回繰り返し、θ̂1∗ , . . . , θ̂Bする。
(iii). 偏り修正推定量 zˆ0 と、加速定数推定量 â を、式 (11)、√(iii). ブートストラップ t 値 T ∗ = n(θ̂∗ − θ̂)/σ̂ ∗ を計(12) から計算する。
∗∗∗∗
算し、B 回繰り返し、T1 , . . . , TB を求める。
(iv). (θ̂α , θ̂1−α ) = (θ̂(B) を θ に対する信頼度αb) , θ̂(B 1[−α)∗∗(iv). wα = T(Bα) とする。Bα は自然数とし、T(Bα) はパー1−2α の両側信頼区間とする。ここで、任意の α(0 ≤センタイル法の時と同様の順序統計量とする。このα ≤ 1) に対して、α̂ は次のように求める。
とき (θ̂ − n−1/2 σ̂w1−α , θ̂ − n−1/2 σ̂wα ) は、名目上のẑ0 + zα
被覆率 1 − 2α をもつ両側信頼区間である。
α̂ = Φ(ẑ0 +)(13)
1 − â(ẑ0 + zα )3.4 BCa 法
1\− α は、式 (13) の zα を z1−α で置き換えたもので推定量 θ̂ は、偏りを持っているかもしれないし、θ̂ の分ある。
布は歪みを持っているかもしれない。パーセントタイル法では θ̂ に対するブートストラップ推定量 θ̂∗ の分布を修 4 データについて正せずにパーセント点を求めているため、推定量の偏り[1] の P84 表 3 のデータから取られる平均と分散共分散や分布の歪みの影響を受けてしまう。そのため、偏りや行列を用いる。
分布の歪みの可能性がある時は、パーセントタイル法で[1] から得られた平均と分散共分散行列を母集団の分散信頼区間を構成するのはあまり良くない。
共分散行列と平均として、mvrnorm でサンプルデータをBCa 法は、推定量の偏りとその分布の歪みを同時に修作成し、各信頼区間を求める。サンプルサイズは、50、正する方法である。この方法を適用するためには、偏り200、500 の 3 種とする。また、平均と分散共分散行列、修正量 z0 と加速定数と呼ばれる量 a を推定しなければな固有値と固有ベクトルは下の表のものである。
らない。
θ̂∗ をパーセンタイル法で定義した量とする。この時、元平均のデータから計算される θ̂ の値は、θ̂∗ の分布の中央に位置
14938.772.23333 79.36667
するものもあれば、そうでないものもある。この食い違い
分散共分散行列を表す量 θ̂∗ の分布の中心と θ̂ との偏差 z0 = Φ−1 (P r[θ̂∗ ≤53.51724 40.79310 27.58621 28.75862θ̂]) は、推定値 θ̂ の θ に対する偏りの程度を表す尺度と考40.79310 41.73448 29.83103 24.35517
えられる。Φ は標準正規分布の分布関数である。ブート
27.58621 29.83103 26.52989 17.22184
ストラップ推定値 θ̂∗ のちょうど半分が θ̂ より小さい時、28.75862 24.35517 17.22184 18.24023
P r[θ̂∗ ≤ θ̂] ＝ 0.5 となり、z0 = 0 となる。z0 は、ブート
ストラップ推定値 θ̂∗ = θ(Y1∗ , . . . , Yn∗ ) を用いて
表 1 母集団の平均と分散共分散行列ẑ0 = Φ−1 (#{θ̂b∗ ≤ θ̂}/B)(11)と表わされる。
この偏り修正量の推定値 ẑ0 が信用できるものであるた
めには、式 (11) における B をパーセンタイル法の B より、倍程度大きくとる必要がある。
加速定数 a は、ブートストラップ標本に基づいて推定を行える。しかし、計算量を減らすために、正則な推定量 θ̂ に対して次のような処理を行う。
θ̂(i) を、元のデータから i 番目のデータ yi を取り除いたものから計算された θ̂ の値とする。ここで、θ̂(・) =∑nn−1 i=1 θ̂(i) と定義すると、加速定数 â は∑n(θ̂(・) − θ̂(i) )3â = ∑ni=1
6{ i=1 (θ̂(・) − θ̂(i) )2 }3/2となる。

124.814224-0.6240226-0.5591912-0.4083340-0.3621661

1.863396-0.37924835-0.10801960-0.08414179
0.91510798

表 2 母集団の固有値と固有ベクトル5(12)固有値10.848728
2.495491
固有ベクトル
0.6455638
0.22363789-0.3456392 -0.74576837-0.6604700 0.62447014
0.1660133
0.06206998

シミュレーションについてシミュレーションによる性能の評価基準として、被覆率を使う。

5.1

被覆率についての被覆率の結果である。

図 2 の縦軸は頻度を、横軸は第 1 固有値、(1,1) 固有ベクトルの大きさを示している。分布の両側にある縦線は両側信頼限界である。母集団の第 1 固有値は 124.814224であり、母集団の (1,1) 固有ベクトル −0.6240226 である。
母集団の固有値、固有ベクトルが両側信頼限界の間に入る時、それぞれ被覆回数を+1 する。被覆率は、図 2 の様なことを 1000 回繰り返して、母集団の値が含まれる割合を調べるものである。

固有値漸近法.945
対数を用いた漸近法 .943

.964.956.969.973.949.940

表 4 漸近法固有値
漸近法固有ベクトル.946 .959 .919 .967.952 .941 .820 .915
.954 .949 .857 .913
.960 .954 .908 .981
表 5 漸近法固有ベクトル

.957
.941
.944
.934.935

固有値.938 .956
固有ベクトル.947 .948.939 .905.949 .917.949 .948.877.952.953.944.945.827.940.943.936.934

表 6 パーセンタイル法

7

図 2 両側信頼区間と第 1 固有値、(1,1) 固有ベクトル
漸近固有値
漸近対数固有値
パーセンタイル法
t法
BCa 法

下側信頼限界112.6412
111.2681103.6404
110.3003
106.2877

上側信頼限界157.1251154.7389
158.4347
174.1243
160.8743

表 3 第 1 固有値両側信頼限界

固有値.940 .912
固有ベクトル.946 .957.938 .885.947 .938.946 .930.903.939.949.937.858

表 7 BCa 法

考察

シミュレーションの結果から、サンプルが正規分布の時、固有値は、n が大きいと全ての手法の差異が小さくなり、固有ベクトルは、特定の要素を除いて、n の大きさにかかわらず全ての手法が似たような結果を出すという事が分かった。

8

おわりに今後の課題としては、サンプルデータを正規分布以外のもので行うことが考えられる。これらを行う事によって、BCa 法と t 法、パーセンタイル法の性能の比較を行い、計算時間を考慮した利便性と正確性の比較も考えられる。

表 3 は、各手法の第 1 固有値の両側信頼点の 1 例を示 参考文献している。サンプルサイズ n は 200 としている。また、 [1] 田中豊、脇本和昌：多変量統計解析法, 現代数学社
信頼区間は 95% 、パーセンタイル法、t 法、BCa 法は、(1983).
ブートストラップ反復回数を 1000 回としている。
[2] 北川源四郎、竹村彰通：21 世紀の統計科学 3:数理・計算の統計科学, 東京大学出版会 (2008).
5.2 プログラムについて
ブートストラップ法と漸近法の固有ベクトルは向きの [3] 伊庭幸人、汪金芳、田栗正章、手塚集、樺島祥介、上田修功：統計科学のフロンティア 11-計算統計 , 岩
制限をかけるため、母集団の固有ベクトルとサンプルの波書店 (2003).
固有ベクトルの内積が正になるように取っている。
t 法ではジャックナイフ分散推定をブートストラップ反 [4] 前園宜彦：統計的推測の漸近理論, 九州大学出版会
復回数分行っており、被覆率を出すのに非常に時間がか(2001).
かるため、t 法のシミュレーションを行わない事にした。 [5] I.T. Jolliﬀe：Principal Component Analysis (Springer
Series in Statistics) , Springer(2002).
6 各手法の被覆率
この章の各表は被覆率を計算する回数を 1000 回、サン
プルサイズ n=200、ブートストラップ反復回数 1000 回[6] Efron.B,Tibshirani,R.J.：An Introduction to theBootstrap. Chapman & Hall,New York(1993).