高次元データの統計的方法論

日本統計学会誌
第 43 巻, 第 1 号, 2013 年 9 月
123 頁 ∼ 150 頁
日本統計学会研究業績賞受賞者特別寄稿論文
高次元データの統計的方法論
青嶋 誠∗ ，矢田 和善†
Effective Methodologies for High-Dimensional Data
Makoto Aoshima∗ and Kazuyoshi Yata†
本論文は，高次元データの統計的方法論を考える．まず，高次元データに対する従来型の主成
分分析（PCA）の限界を示す．それに替わるものとして「ノイズ掃き出し法」と「クロスデータ
行列法」という 2 つの新しい PCA 手法を紹介する．高次元の推測に現れる各種パラメータに対
して，
「拡張クロスデータ行列法」という推定法を紹介し，これが計算コストを著しく削減し，且
つ，漸近分散が小さな不偏推定量を与えることを示す．さらに，高次元における相関ベクトルや
母集団が複数ある場合の平均ベクトルに関する推定と検定を紹介し，推測において高い精度を保
証するための標本数の決定法を与える．最後に，高次元データの判別分析について，母集団間の
ユークリッド距離に基づく判別方式と幾何学的差異を利用した判別方式を紹介し，それらが誤判
別確率に関して高い精度を保証することを示す．
In this paper, we consider statistical methodologies for high-dimensional data. We first
clarify the limit of the conventional principal component analysis (PCA) for high-dimensional
data. In order to overcome the curse of dimensionality, we introduce two effective PCAs called
the noise-reduction methodology and the cross-data-matrix (CDM) methodology. We further
introduce the extended CDM methodology, which offers an unbiased estimator having small
asymptotic variance and low computational cost, for parameters appearing in high-dimensional
data analysis. We give correlation tests and inference on multiclass mean vectors for highdimensional data, and discuss sample size detemination to ensure prespecified high accuracy for
inference. Finally, we provide two effective discriminant procedures: a distance-based classifier
and a geometric classifier, which can ensure high accuracy in misclassification rates.
キーワード: 幾何学的表現，クラスター分析，高次元小標本，主成分分析，2 段階検定，パス解
析，判別分析，標本数決定
はじめに
1.
ゲノム科学，情報工学，金融工学などの現代科学の 1 つの特徴は，データがもつ次元数
の膨大さにある．例えば，次世代シークエンサによるゲノム配列データなど，次元数が数
百万を超えるデータも解析の対象になる．こういった高次元データの第 1 の特徴は，次元
∗
†
筑波大学数理物質系：〒 305-8571 茨城県つくば市天王台 1-1-1 (E-mail: [email protected])．
筑波大学数理物質系：〒 305-8571 茨城県つくば市天王台 1-1-1 (E-mail: [email protected])．
124
日本統計学会誌 第43巻 第1号 2013
数が標本数を遥かに超えることである．第 2 の特徴は，高次元データは豊富な情報を有す
るものの，それが巨大なノイズに埋もれ見つけ難いことである．これらの理由から，通常
の多変量解析法では高次元データの推測に精度を保証することができず，間違った解析結
果を導くことさえある．そのため，高次元データの解析には，新しい理論と方法論が必要
になる．大標本漸近理論だけでは十分とはいえず，母集団に正規分布を仮定した方法論で
は厳しいといわざるを得ない．高次元データの解析には，高次元漸近理論の構築が必要に
なり，母集団に正規分布を仮定しない方法論，もしくは，ノンパラメトリックな方法論が
必要になる．さらに，膨大なデータを処理するためには，低い計算コストで高精度な解析
結果を出力できるようなアルゴリズムが求められる．
上で述べたような理論と方法論を構築するためには，高次元データのもう 1 つの特徴に
目を向けることが重要になる．それは，データの次元数が標本数を遥かに超えた高次元空
間に現れるデータの幾何学的特徴である．Ahn et al. (2007)，Hall et al. (2005) は，標本
数 n を固定して次元数 p を p → ∞ としたときの高次元小標本について，高次元データの
幾何学的表現を研究した．ただし，彼らの研究は，母集団に正規分布もしくは類する条件
を仮定する限定的なものであった．Yata and Aoshima (2012a) は，高次元小標本の母集団
が非正規分布の場合を研究し，先行研究では見つけることができなかった高次元データの
2 つの幾何学的表現を発見した．非正規性の尺度を境にした，固有空間の球面集中現象と
座標軸集中現象である．Aoshima and Yata (2011a) は，高次元データの統計的推測に，幾
何学的表現に基づいた各種方法論を考案し，統計量の高次元漸近正規性，標本数の設計，推
測の精度保証に至るまで，一連の基礎理論を築いた．Aoshima and Yata (2011a) の研究は
多岐に亘り，与えられたバンド幅をもつ高次元球面上の信頼領域，高次元二標本問題，高
次元共分散行列の推定・検定，高次元判別分析，高次元回帰分析，高次元変数選択，パス
ウエイ解析など，様々な統計的推測の問題に対して，高次元データの豊富な情報を生かす
ための理論と方法論を考えている．本稿では，Aoshima and Yata の一連の研究で得られ
た高次元データの方法論について，その有用性と使い方を解説する．高次元データの理論
に関する発展については，青嶋・矢田 (2013) を参照されたい．
まず，第 2 節で，Yata and Aoshima (2009) に基づいて，高次元データに対する従来
型の主成分分析 (PCA) の使用限界を解説する．その際に，確率論と理論物理の方面から
Johnstone (2001) 等によって提案された固有値モデルが，実際の高次元データと乖離する
ことを示す．第 3 節では，Yata and Aoshima (2010a, b, 2012a) が高次元データに対して
考案した「ノイズ掃き出し法」と「クロスデータ行列法」という 2 つの新しい PCA を紹介
する．その際に，これら方法論のアイディアの基となる高次元データの 2 つの幾何学的表
現を示す．第 4 節では，Aoshima and Yata (2013a)，Yata and Aoshima (2013a) に基づい
て，高次元の推測に現れる各種パラメータに対して「拡張クロスデータ行列法」という推定
125
高次元データの統計的方法論
法を紹介する．これが計算コストを著しく削減し，且つ，漸近分散が小さな不偏推定量を
与えることを示す．さらに，高次元における相関ベクトルの検定を考え，拡張クロスデー
タ行列法で構築した統計量が高次元漸近正規性を有すること，及び，推測に高い精度を保
証するための標本数決定法を示す．第 5 節では，Aoshima and Yata (2011a, b, 2013b) に
基づき，母集団が複数ある場合の平均ベクトルに関して高次元における推定と検定を論じ
る．最後に，第 6 節では，高次元データの判別分析について，Aoshima and Yata (2011a,
2013a) によって提案された，母集団間のユークリッド距離に基づく判別方式と幾何学的差
異を利用した判別方式を紹介する．それらが誤判別確率に関して高い精度を保証すること
を示す．なお，本論文で紹介する方法論は，母集団に正規分布の仮定を必要としない．し
かも，高次元データの様々な推測に高精度を保証する，簡便かつ高速で頑健な方法論であ
る．高次元データのみならず，膨大で多様なデータを含むビックデータにも対応し得る，
汎用性を有した方法論である．
2.
高次元データに対する従来型 PCA の限界
平均に p 次のベクトル µ(p) ，共分散行列に p 次の正定値対称行列 Σ(p) (> O) をもつ母集
団を考える．母集団から n (≥ 4) 個の p 次データベクトル x1(p) , . . . , xn(p) を無作為に抽出
し，p × n データ行列 X (p) = [x1(p) , . . . , xn(p) ] を定義する．Σ(p) の固有値を λ1(p) ≥ · · · ≥
λp(p) (> 0) とし，適当な直交行列 H (p) = [h1(p) , . . . , hp(p) ] で Σ(p) を Σ(p) = H (p) Λ(p) H T(p) ,
−1/2
Λ(p) = diag(λ1(p) , . . . , λp(p) ) と分解する．そのとき Z (p) = Λ(p) H T(p) (X (p) − [µ, . . . , µ])
とおき，Z (p) = [z 1(p) , . . . , z p(p) ]T ，z i(p) = (zi1(p) , . . . , zin(p) )T と表記する．ただし，Z (p)
の成分は，4 次モーメントに一様有界性を仮定する．今後，次元数を意識して付した添え
字 (p) は省略する．本節では，必要であれば Z に以下を仮定する．
(A-i) zij , i = 1, . . . , p が互いに独立．
(A-i) は，X の母集団分布について，正規分布よりも緩い仮定になっている．
高次元データ空間の特徴を捉えるために，有名な 3 つの DNA マイクロアレイデータを紹
介する．1 つ目は，Alon et al. (1999) の結腸データで，遺伝子数は 2000，腫瘍患者 (Colon)
の 40 サンプルと正常者 (Normal) の 22 サンプルからなる．2 つ目は，Golub et al. (1999)
の白血病データで，遺伝子数は 7129，ALL タイプの 47 サンプルと AML タイプの 25 サン
プルからなる．3 つ目は，Singh et al. (2002) の前立腺データで，遺伝子数は 12600，正常
者 (Normal) の 50 サンプルと腫瘍患者 (Prostate) の 52 サンプルからなる．各々のデータ
セットを標準化し，第 3 固有値まで推定したものが表 1 である．λ̂i は従来型の方法による
推定値を表し，後述する (2.3) 式から計算した．一方，λ́i と λ̃i は，第 3 節で紹介するノイ
ズ掃き出し法とクロスデータ行列法による推定値を表し，それぞれ (3.7) 式と (3.8) 式から
126
表1
日本統計学会誌 第43巻 第1号 2013
3 つの DNA マイクロアレイデータの第 3 固有値までの推定値．λ̂i は従来型の推定値，λ́i はノイズ掃き
出し法による推定値，λ̃i はクロスデータ行列法による推定値を表す．
n
n/p
λ̂1 , λ̂2 , λ̂3
λ́1 , λ́2 , λ́3
λ̃1 , λ̃2 , λ̃3
遺伝子数 2000 (= p) の結腸データ (Alon et al. (1999))
Colon
40
0.02
949, 228, 150
922, 205, 131
895, 194, 98
Normal
22
0.011
922, 170, 137
868, 122, 94
827, 137, 91
遺伝子数 7129 (= p) の白血病データ (Golub et al. (1999))
ALL
47
0.0066
1148, 841, 342
1015, 724, 231
941, 685, 208
AML
25
0.0035
1344, 733, 488
1093, 504, 271
1004, 441, 260
遺伝子数 12600 (= p) の前立腺データ (Singh et al. (2002))
図1
Normal
50
0.004
6748, 561, 426
6626, 448, 320
6360, 331, 287
Prostate
52
0.0041
6095, 685, 511
5965, 566, 401
5987, 568, 370
表 1 の 6 つのデータセットに対する，ノイズ掃き出し法による第 15 固有値までの推定値．
計算した．図 1 では，ノイズ掃き出し法による第 15 固有値までの推定値をプロットした．
表 1 と図 1 から，各データは標本数と次元数の比 n/p が非常に小さいこと，そして，最初
のいくつかの固有値が飛び抜けて大きく，他の固有値は相対的に小さなものであることが
分かる．さらに，データセット間で比較すると，最初のいくつかの固有値は次元数に関し
て指数関数的に大きなものであることが伺える．
Σ の固有値について，Johnstone (2001) は次のような spiked モデルを仮定した．
λi > 1 (i = 1, . . . , m) は p に依存しない定数,
λi = 1 (i = m + 1, . . . , p).
(2.1)
このモデルのもと，p と n が同等なオーダーで発散する高次元大標本 (n/p → c > 0) の
枠組みで，多くの研究者が固有値の漸近的な振る舞いを研究をした．例えば，Johnstone
高次元データの統計的方法論
127
(2001)，Baik et al. (2005)，Paul (2007)，Johnstone and Lu (2009) は母集団に正規分布を
仮定し，また，Baik and Silverstein (2006)，Lee et al. (2010) は前述の (A-i) よりも限定
的な条件を仮定して，固有値の漸近理論を研究している．固有値モデル (2.1) の後半 (λi の
m + 1 番目以降) はノイズを表す部分である．しかし，すべてのノイズが等しいという設定
は，数学的な扱いを簡単にするものの，いささか厳しい制約である．一方で，(2.1) の前半
は潜在的な固有空間を表す部分である．表 1 と図 1 を見ると，固有値が次元数 p に依存しな
いという設定は，実際のデータから乖離していると言わざるを得ない．Yata and Aoshima
(2009) は，次のような一般化 spiked モデルを考えた．
λi = ai pαi (i = 1, . . . , m),
λi = ci (i = m + 1, . . . , p).
(2.2)
ここで，ai (> 0), ci (> 0), αi (α1 ≥ · · · ≥ αm > 0) はともに p に依存しない未知の実数，
m は p に依存しない未知の自然数で，これらのパラメータは λ1 ≥ · · · ≥ λp なる大小関係
を満たすものとする．
注意 1
一般化 spiked モデル (2.2) は，高次元データ空間の固有値の大きさと，それを推
定するための標本数と次元数の関係を見るための，簡便なモデルである．Yata and Aoshima
(2013b) は，(2.2) を拡張したモデルとして power spiked モデルを提案し，高次元 PCA の
漸近理論を展開している．
∑n
標本共分散行列を S n = (n−1)−1 (X −X)(X −X)T = (n−1)−1 i=1 (xi − x̄)(xi − x̄)T
∑n
とする．ここで，X = [x̄, . . . , x̄]，x̄ = j=1 xj /n である．そのとき，n × n の行列 S nD =
(n − 1)−1 (X − X)T (X − X) を S n と正の固有値を共有する双対標本共分散行列という．
S nD の固有値を λ̂1 ≥ · · · ≥ λ̂n−1 (≥ 0)，λ̂j に対する固有ベクトルを ûj とし，スペクトル
分解
S nD =
n−1
∑
λ̂j ûj ûTj
(2.3)
j=1
を使って従来型 PCA を定義する．Jung and Marron (2009) は，母集団に正規分布もしく
は ρ-mixing を仮定して，p → ∞ で n は固定のもと，(2.2) における αi (> 1) をもつ λi に
対して推定量 λ̂i の一致性を論じた．ところが，αi > 1 なる状況は，数学的な扱いは簡単で
あるものの，実際のデータ解析ではあまり意味がない．Jung et al. (2012) は，p → ∞ で n
は固定のもと，αi = 1 の場合を扱った．一方，Yata and Aoshima (2009) は，母集団の分
布型や αi に制約を課すことなしに，λ̂i が一致性をもつための条件を導いた．標本数 n が
次元数 p に依存するとき n(p) で表し，n(p) = pγ (γ > 0) とおく．そのとき，例えば固有
値について，Yata and Aoshima (2009) は次の 2 つの定理を与えた．
128
図2
日本統計学会誌 第43巻 第1号 2013
仮定 (A-i) とモデル (2.2) のもと，従来型の推定量 λ̂i が一致性を有する領域 (白：γ > 1 − αi ) と一致性
を有さない領域 (灰色：γ ≤ 1 − αi )．横軸は Spike Index αi (λi = ai pαi )，縦軸は Index γ (n = pγ )．
定理 2.1
(A-i) を仮定する．(2.2) に関して条件
(i) αi > 1 のとき p → ∞, n → ∞,
(ii) αi ∈ (0, 1] のとき p → ∞, p1−αi /n(p) → 0
を満たせば，λ̂i (i = 1, . . . , m) は次が成り立つ．
λ̂i
= 1 + oP (1).
λi
定理 2.2
(2.4)
(A-i) を仮定しない．(2.2) に関して次の条件を満たせば，(2.4) が成り立つ．
(i) αi > 1 のとき p → ∞, n → ∞,
(ii) αi ∈ (0, 1] のとき p → ∞, p2−2αi /n(p) → 0.
これらの定理は，従来型の固有値推定が一致性をもつための条件を与える．固有値の冪が
αi ∈ (0, 1] のとき，標本数 n は次元数 p に依存して決めるべきだと分かる．図 2 は，(A-i)
を仮定するときに，Yata and Aoshima (2009) で一致性が証明された領域 αi < 1 と，Jung
and Marron (2009) と Jung et al. (2012) で一致性が証明された領域 αi ≥ 1 を示したもの
である．例えば，マイクロアレイデータを想定して，p = 10000, n = 100 (= p1/2 ) とする．
そのとき，従来型の推定量が一致性を有する領域は，(A-i) を仮定する場合は定理 2.1(ii) か
ら αi > 1/2，(A-i) を仮定しない場合には定理 2.2(ii) から αi > 3/4 となって，推定対象と
なり得る固有値は厳しく限定される．なお，(A-i) を仮定しない例は 3.1 節の図 4 を参照の
こと．Yata and Aoshima (2009) は，固有ベクトルや主成分スコアの従来型推定について
も一致性の条件を与えている．次節では，高次元データに対する PCA の限界を緩めるた
めに，Yata and Aoshima (2010a, b, 2012a) が考案した 2 つの方法論を紹介する．
129
高次元データの統計的方法論
ノイズ掃き出し法とクロスデータ行列法
3.
第 2 節と同様に，平均に p 次のベクトル µ，共分散行列に p 次の正定値対称行列 Σ (> O)
をもつ母集団を考え，Z の成分に 4 次モーメントの一様有界性を仮定する．高次元小標本
の幾何学的表現に基づいて，
「ノイズ掃き出し法」と「クロスデータ行列法」という 2 つの
PCA を紹介する．
3.1
高次元小標本の幾何学的表現
高次元データを解析する上で鍵となるのは，高次元小標本の幾何学的表現である．表記
を簡単にするために，µ = 0 とする．そのとき，双対標本共分散行列を S onD = n−1 X T X
∑n
とおき，スペクトル分解を S onD = j=1 λ̂oj ûoj ûToj (λ̂o1 ≥ · · · ≥ λ̂on ) とする．固有値に
次のような球形条件を考える．
∑p
λ2
∑pi=1 i 2 → 0,
( i=1 λi )
p → ∞.
(3.1)
そのとき，Yata and Aoshima (2012a) は，高次元小標本の幾何学的表現について次の定理
を与えた．
定理 3.1
球形条件 (3.1) を仮定する．もしも，Z が
∑p
2
2
i,j≥1 λi λj E{(zik − 1)(zjk − 1)}
∑p
→ 0,
( j=1 λj )2
p→∞
(3.2)
なる条件を満たせば，p → ∞（n は固定したままでよい）で次が成り立つ．
n
∑p
P
i=1 λi
S onD −
→ I n;
(3.3)
wj ∈ {en ∈ Rn | k en k= 1},
ここで，I n は n 次単位行列，wj = (n/
∑p
i=1
j = 1, . . . , n.
λi )λ̂oj ûoj である．一方，Z が (3.2) を満た
さない場合，p → ∞（n は固定したままでよい）である正則条件のもと次が成り立つ．
n
∑p
i=1 λi
P
S onD −
→ Dn ;
ûoj ∈ {e1n , · · · , enn },
(3.4)
j = 1, . . . , n.
ここで，D n は対角成分が OP (1) となる対角行列，ein (i = 1, . . . , n) は第 i 成分が 1 で他
は 0 となる n 次基本ベクトルである．
幾何学的表現の理論的背景と意味については，青嶋・矢田 (2013) を参照されたい．(3.2) を
満たす例と満たさない例として，p 次元正規分布 Np (0, I p ) と，平均が 0 で共分散行列が I p で
自由度が 5 の p 次元 t 分布を考える．標本数は n = 3 とする．次元数が p = 4, 40, 400, 4000
130
日本統計学会誌 第43巻 第1号 2013
図3
(a) p = 4
(b) p = 40
(c) p = 400
(d) p = 4000
p 次元正規分布 Np (0, I p ) における 15 組の ±wj (j = 1, 2, 3) の幾何学的表現．
の各場合で，±wj (j = 1, 2, 3) の対を独立に 15 組発生させて，出力結果を w1 は
は
で，w3 は
で，w2
で表示する．図 3 の (a)–(d) は正規分布について，図 4 の (a)–(d) は t
分布について，結果を纏めたものである．次元数が増えるにつれて，(3.2) を満たす場合と
そうでない場合の幾何学的差異が，明瞭になっていく様子が見てとれる．高次元小標本が
これらの幾何学的表現を有することを承知しておくことは，高次元データに対する理論と
方法論を構築する上で重要といえる．例えば，S onD は
nS onD =
m
∑
j=1
λj z j z Tj
+
p
∑
λj z j z Tj
j=m+1
と分解できる．この分解の第 2 項は，固有値モデル (2.2) のもとでノイズを表し，球形条件
(3.1) を満たす．母集団に (A-i) を仮定する．その場合，定理 3.1 の (3.3) 式より
∑p
T
j=m+1 λj z j z j P
∑p
−
→ I n, p → ∞
j=m+1 λj
(3.5)
131
高次元データの統計的方法論
図4
(a) p = 4
(b) p = 40
(c) p = 400
(d) p = 4000
平均 0，共分散行列 I p ，自由度 5 の p 次元 t 分布における 15 組の ±wj (j = 1, 2, 3) の幾何学的表現．
なる幾何学的表現を得る．(3.5) 式を S onD の分解に適用すれば，
ûToi S onD ûoi


p
m
T
∑
∑
λ
z
z
λj
j j j
 ûoi +
= ûToi 
{1 + op (1)},
n
n
j=m+1
j=1
p→∞
(3.6)
なる分解が得られ，固有値の推定量に含まれるノイズの主要項が見積もれる．
3.2
ノイズ掃き出し法による新しい PCA
母集団に (A-i) を仮定できる場合に，Yata and Aoshima (2012a) はノイズ掃き出し法と
いう方法論を提案した．(A-i) を仮定すれば，図 3 に見られるような幾何学的表現をもち，
固有値モデル (2.2) においては (3.5) 式が成り立つ．Yata and Aoshima (2012a) は，(3.6)
式のノイズの主要項を掃き出すために，次のような固有値の推定量を提案した．
λ́i = λ̂i −
tr(S nD ) −
∑i
j=1
n−1−i
λ̂j
(i = 1, . . . , n − 2).
(3.7)
132
図5
日本統計学会誌 第43巻 第1号 2013
仮定 (A-i) とモデル (2.2) のもと，ノイズ掃き出し法による推定量 λ́i が一致性を有する領域 (γ > 1 − 2αi )
と一致性を有さない領域 (灰色：γ ≤ 1 − 2αi )．横軸は Spike Index αi (λi = ai pαi )，縦軸は Index γ (n = pγ )．
(3.7) 式の第 2 項でノイズの主要項が掃き出される．これが，ノイズ掃き出し法である．Yata
and Aoshima (2012a) は次の定理を与えた．
定理 3.2
(A-i) を仮定する．(2.2) に関して条件
(i) αi > 1/2 のとき p → ∞, n → ∞,
(ii) αi ∈ (0, 1/2] のとき p → ∞, p1−2αi /n(p) → 0
を満たせば，λ́i (i = 1, . . . , m) は次が成り立つ．
λ́i
= 1 + oP (1).
λi
図 5 は，ノイズの掃き出し法による推定量が一致性を有する領域を示す．図 2 と比べる
と，従来型 PCA よりも一致性を有する領域が広がっていることが分かる．ノイズの主要
∑p
項 j=m+1 λj /n を掃き出したことで，一致性をもつための標本数 n の条件が軽減された
のである．Yata and Aoshima (2012a) は次の定理を与えた．
定理 3.3
(A-i) を仮定する．(2.2) における単根の固有値 λi (i = 1, . . . , m) に関して
条件
(i) αi > 1/2 のとき p → ∞, n → ∞,
(ii) αi ∈ (0, 1/2] のとき p → ∞, p2−4αi /n(p) → 0
を満たせば，λ́i (i = 1, . . . , m) は次が成り立つ．
(
)
√
n
λ́i
− 1 ⇒ N (0, 1).
M i λi
133
高次元データの統計的方法論
2
ここで，Mi = Var(zij
) (> 0), i = 1, . . . , p (j = 1, . . . , n)，⇒ は分布収束を表す．
次に，ノイズ掃き出し法による主成分スコアの推定を考える．データ xj の第 i 主成分ス
√
コアは，hTi (xj − µ) = λi zij (= sij とおく) である．(2.3) 式における S nD の固有ベク
トルの成分を ûi = (ûi1 , . . . , ûin )T とする．推定量 (3.7) に基づいて，第 i 主成分スコアを
√
ûij nλ́i (= śij とおく) で推定する．Yata and Aoshima (2012a) は次の定理を与えた．
定理 3.4
(A-i) を仮定する．MSE(śi ) = n−1
∑n
j=1 (śij
− sij )2 とおく．(2.2) における
単根の固有値 λi (i = 1, . . . , m) に関して定理 3.2 の条件 (i)–(ii) を満たせば，MSE(śi ) につ
いて次が成り立つ．
MSE(śi )
= oP (1).
λi
(A-i) が仮定されるとき，ノイズ掃き出し法による固有値，固有ベクトル，主成分スコア
の推定は，従来型 PCA を著しく改良することが報告されている．
3.3
クロスデータ行列法による新しい PCA
母集団に (A-i) が仮定できない場合，Yata and Aoshima (2010a) は，母集団分布の仮定
を必要としないクロスデータ行列法という方法論を考案した．この場合，図 4 のような座
標軸集中現象も考慮に入れなければならず，ノイズの大きさが定まらないために，問題は
格段に難しくなる．
クロスデータ行列法は次のような方法論である．データ行列 X を（無作為に）2 つに分割
して，p × n(l) 部分行列 X l = [xl1 , · · · , xln(l) ] (l = 1, 2) を定義する．ここで，n(1) = dn/2e，
∑n(l)
n(2) = n − n(1) ，dxe は x 以上の最小の整数を表す．X l = [x̄l , . . . , x̄l ], x̄l = n−1
j=1 xlj
(l)
(l = 1, 2) とおく．そのとき，クロスデータ行列を S nD(1) = {(n(1) − 1)(n(2) − 1)}−1/2 (X 1 −
X 1 )T (X 2 − X 2 ), S nD(2) = S TnD(1) と定義し，S nD(1) の特異値分解を
n(2) −1
S nD(1) =
∑
λ̃i ũi(1) ũTi(2)
(3.8)
i=1
とおく．ここで，λ̃1 ≥ · · · ≥ λ̃n(2) −1 (≥ 0) は S nD(1) の特異値，ũi(1) は左特異ベクトル，
ũi(2) は右特異ベクトルである．特異値分解 (3.8) に基づく方法論をクロスデータ行列法と
よぶ．例えば，Σ の固有値 λi と固有ベクトル hi は，次のように推定する．
[クロスデータ行列法による固有値・固有ベクトルの推定]
(手順 1) X = [x1 , . . . , xn ] を無作為に 2 分割し，X 1 = [x11 , . . . , x1n(1) ], X 2 = [x21 , . . . ,
x2n(2) ] とおく．
134
日本統計学会誌 第43巻 第1号 2013
(手順 2) クロスデータ行列 S nD(1) = {(n(1) − 1)(n(2) − 1)}−1/2 (X 1 − X 1 )T (X 2 − X 2 ),
S nD(2) = S TnD(1) を定義する．
(手順 3) 対称行列 S nD(1) S nD(2) の固有値 λ̃21 ≥ · · · ≥ λ̃2n(2) −1 (≥ 0) と，対応する固有ベ
√
クトル ũi(1) (i = 1, . . . , n(2) − 1) を計算する．固有値 λi を特異値 λ̃i = λ̃2i で推定
する．
(手順 4) 手順 3 と同様に，対称行列 S nD(2) S nD(1) の固有ベクトル ũi(2) (i = 1, . . . , n(2) −1)
を計算し，ũi(2) = Sign(ũTi(1) S nD(1) ũi(2) )ũi(2) で符号を調整する．
−1/2
(手順 5) S nD(1) の特異値 λ̃i と特異ベクトル ũi(l) , l = 1, 2 に基づいて h̃i = λ̃i
{(n(1) −
1)−1/2 (X 1 − X 1 )ũi(1) + (n(2) − 1)−1/2 (X 2 − X 2 )ũi(2) }/2 を計算し，固有ベクトル hi
を h̃i∗ = h̃i / k h̃i k で推定する．
注意 2
X 1 と X 2 の分割は n Cn(1) 通りある．理論上は無作為に分割してよいが，標本
数 n が非常に小さいと，結果が分割に影響される場合もある．Aoshima and Yata (2011c)
は，2 分割を計算機上でランダムに実行する「一般化クロスデータ行列法」を与えた．
クロスデータ行列 S nD(1) の特異値分解に着目する理由は，高次元データの次のような
幾何学的特徴にある．X の分割に対応して Z を 2 つに分割し，z 1j = (z1j1 , . . . , z1jn(1) )T ，
∑n(l)
z 2j = (z2j1 , . . . , z2jn(2) )T , j = 1, . . . , p を定義する．そのとき，z̄lj = n−1
k=1 zljk , (j =
(l)
1, . . . , p; l = 1, 2) とし，z olj = z lj − (z̄lj , . . . , z̄lj )T (j = 1, . . . , p; l = 1, 2) とおくと，
{(n(1) − 1)(n(2) − 1)}1/2 S nD(1) =
m
∑
λj z o1j z To2j +
j=1
p
∑
λj z o1j z To2j
j=m+1
と表せる．この第 2 項は，固有値モデル (2.2) のもとでノイズを表し，無条件で
∑p
T
j=m+1 λj z o1j z o2j
∑p
j=m+1 λj
P
−
→ O,
p→∞
なる幾何学的特徴をもつ．これは，クロスデータ行列を用いれば，ノイズを自動的に除去
できることを意味する．Yata and Aoshima (2010a) は，ノイズの漸近的挙動を精密に評価
して次の 2 つの定理を与えた．
定理 3.5
(2.2) に関して条件
(i) αi > 1/2 のとき p → ∞, n → ∞,
(ii) αi ∈ (0, 1/2] のとき p → ∞, p2−2αi /n(p) → 0
高次元データの統計的方法論
135
を満たせば，λ̃i (i = 1, . . . , m) は次が成り立つ．
λ̃i
= 1 + oP (1).
λi
定理 3.6
(2.2) における単根の固有値 λi (i = 1, . . . , m) に関して定理 3.5 の条件 (i)–(ii)
を満たせば，λ̃i (i = 1, . . . , m) は次が成り立つ．
(
)
√
n
λ̃i
− 1 ⇒ N (0, 1).
M i λi
2
ここで，Mi = Var(zij
) (> 0), i = 1, . . . , p (j = 1, . . . , n) である．
定理 3.5 の条件を定理 2.2 と比べると，λ̃i は従来型の固有値推定よりも緩い条件で一致
性をもつことが分かる．また，定理 3.3 と定理 3.6 を見ると，推定量の漸近正規性が同じ形
で与えられている．つまり，ノイズ掃き出し法による推定量 λ́i とクロスデータ行列法によ
る推定量 λ̃i は漸近的に同等な分散をもつ．Yata and Aoshima (2010a) は，固有ベクトル
の推定量 h̃i∗ についても，従来型の推定量 ĥi = {(n − 1)λ̂i }−1/2 (X − X)ûi よりも緩い条
件で一致性をもつことを証明している．
次に，クロスデータ行列法による主成分スコアの推定を考える．S nD(1) の特異ベクトル
を ũi(l) = (ũi1(l) , . . . , ũin(l) (l) )T (l = 1, 2) と成分表示する．S nD(1) の特異値と特異ベクトル
√
に基づいて，xlj (l = 1, 2) の第 i 主成分スコアを ũij(l) n(l) λ̃i (= s̃ij(l) とおく) で推定す
る．各 i で，s̃i1(1) , . . . , s̃in(1) (1) , s̃i1(2) , . . . , s̃in(2) (2) に s̃ij , j = 1, . . . , n という通し番号を付
ける．Yata and Aoshima (2010a) は次の定理を与えた．
定理 3.7
MSE(s̃i ) = n−1
∑n
j=1 (s̃ij
− sij )2 とおく．(2.2) における単根の固有値 λi (i =
1, . . . , m) に関して定理 3.5 の条件 (i)–(ii) を満たせば，MSE(s̃i ) について次が成り立つ．
MSE(s̃i )
= oP (1).
λi
クロスデータ行列法による固有値，固有ベクトル，主成分スコアの推定は，(A-i) を仮
定できないノンパラメトリックな状況において，従来型の PCA を著しく改良することが
報告されている．クロスデータ行列法の応用は広く，例えば，Yata and Aoshima (2010b)
は高次元データの潜在空間の次元推定に， また，Aoshima and Yata (2011c) と Yata and
Aoshima (2010a) はマイクロアレイデータのクラスター分析に応用している．
4.
拡張クロスデータ行列法
平均に p 次のベクトル µ，共分散行列に p 次の非負定値対称行列 Σ (≥ O) をもつ母集団
を考える．Yata and Aoshima (2013a) は，クロスデータ行列法における 2 分割の組み合わ
136
日本統計学会誌 第43巻 第1号 2013
せを考慮した「拡張クロスデータ行列法」を開発した．本節では，高次元データ解析の鍵
を握る重要なパラメータについて，拡張クロスデータ行列法を用いることで，推定量と検
定統計量が効率よく構築されることを示す．高次元データに次のモデルを考える．
xj = Γwj + µ,
j = 1, . . . , n.
ここで，Γ は ΓΓT = Σ なる p × r 行列，wj は E(wj ) = 0, Var(wj ) = I r とする．このモ
デルは，Bai and Saranadasa (1996)，Chen and Qin (2010)，Aoshima and Yata (2013a)
等で解析された．いま，wj = (w1j , . . . , wrj )T とおき，各成分は 4 次モーメントが一様有
界と仮定する．母集団分布には，必要な箇所でその都度，次のどちらかを仮定する．
(A-ii)
(A-iii)
wij ，i = 1, . . . , r は互いに独立である．
2
2
E(wij
wsj
) = 1，E(wij wsj wtj wuj ) = 0, i 6= s, t, u.
(A-iii) は (A-ii) を緩め，(A-ii) は正規分布を緩めた仮定になっている．
4.1
tr(Σ2 ) の推定量
高次元データの推測に精度を保証するための，鍵となるパラメータの 1 つが tr(Σ2 ) であ
る．単純な推定量 tr(S 2n ) は，高次元データに対して非常に大きなバイアスをもち役に立た
ない．Yata (2010) は，クロスデータ行列法を用いて tr(Σ2 ) の推定を考えた．標本を 2 分割
し，X = [X 1 , X 2 ], X 1 = [x1 , . . . , xn(1) ], X 2 = [xn(1) +1 , . . . , xn ] とする．各分割から標本
共分散行列 S n(l) = (n(l) − 1)−1 (X l − X l )(X l − X l )T , l = 1, 2 を計算する．このとき，不
偏性 E{tr(S n(1) S n(2) )} = tr(Σ2 ) が示され，tr(S n(1) S n(2) ) は tr(Σ2 ) の不偏推定量になる．
母集団分布に (A-iii) を仮定すると，p → ∞, n → ∞ のとき Var{tr(S n(1) S n(2) )/tr(Σ2 )} =
8{1+o(1)}/n2 +O[tr(Σ4 )/{ntr(Σ2 )2 }] → 0 となり，n/p → 0 なる高次元小標本の枠組みで一
致性をもつ．分割 X = [X 1 , X 2 ] に対応したクロスデータ行列 S nD(1) が tr(S nD(1) S TnD(1) ) =
tr(S n(1) S n(2) ) となることに注意する．
c2 ) = c−1 {tr(S 2 )−
一方で，Bai and Saranadasa (1996)，Srivastava (2005) は，推定量 tr(Σ
n
n
tr(S n )2 /(n − 1)} を与えた．ここで，cn = (n − 2)(n + 1)/(n − 1)2 である．そのとき，母
c2 )} = tr(Σ2 ) なる不偏性をもち，p → ∞，n → ∞
集団に正規分布を仮定すると，E{tr(Σ
c2 )/tr(Σ2 )} = 4{1 + o(1)}/n2 + 8tr(Σ4 ){1 + o(1)}/{ntr(Σ2 )2 } → 0 とな
のとき Var{tr(Σ
c2 ) の不偏性は主張できず，高次元
る．しかし，母集団に正規分布を仮定できないと，tr(Σ
において非常に大きなバイアスが生じる．さらに，wj の成分に 8 次モーメントの一様有界
c2 )/tr(Σ2 )} < ∞ さえも保証できない．
性が仮定できないと，Var{tr(Σ
Yata and Aoshima (2013a) は，拡張クロスデータ行列法という方法論を開発した．2 つ
高次元データの統計的方法論
図6
137
b(i + j)/2c ≥ n(1) (i < j) の場合における，集合 V n(1)(i+j) と V n(2)(i+j) のイメージ図．
の集合 V n(1)(k) , V n(2)(k) (k = 3, . . . , 2n − 1) を次のように定義する．


{bk/2c − n(1) + 1, . . . , bk/2c},
bk/2c ≥ n(1) のとき,
V n(1)(k) =

{1, . . . , bk/2c} ∪ {bk/2c + n(2) + 1, . . . , n}, それ以外．


{bk/2c + 1, . . . , bk/2c + n(2) },
bk/2c ≤ n(1) のとき,
V n(2)(k) =

{1, . . . , bk/2c − n(1) } ∪ {bk/2c + 1, . . . , n}, それ以外．
ここで，bxc は x 以下の最大の整数を表す．そのとき，k = 3, . . . , 2n − 1 について，
#(V n(l)(k) ) = n(l) ，l = 1, 2，V n(1)(k) ∩ V n(2)(k) = ∅，V n(1)(k) ∪ V n(2)(k) = {1, . . . , n}
となること，及び，i < j (≤ n) について
i ∈ V n(1)(i+j) ,
j ∈ V n(2)(i+j)
となることに注意する（図 6 参照）
．ここで，#(S) は集合 S の要素の個数を表す．V n(1)(i+j)
と V n(2)(i+j) に対応してデータ集合を 2 分割し，それらに基づく 1 つの不偏推定量を計算し
て，i < j (≤ n) のすべての組合せで平均をとる．これが拡張クロスデータ行列法である．
注意 3
V n(1)(i+j) , V n(2)(i+j) は，計算機の If 関数と Floor 関数で容易に組める．
拡張クロスデータ行列法を使えば，tr(Σ2 ) の不偏推定量は次のように構築できる．各
k (= 3, . . . , 2n − 1) について，2 分割した集合の標本平均を
∑
∑
xn(1)(k) = n−1
xn(2)(k) = n−1
xj ,
xj
(1)
(2)
j∈V n(1)(k)
j∈V n(2)(k)
(4.1)
とし，tr(Σ2 ) の 1 つの不偏推定量として un {(xi − xn(1)(i+j) )T (xj − xn(2)(i+j) )}2 (i < j) を
計算する．ただし，un = n(1) n(2) /{(n(1) − 1)(n(2) − 1)} である．すべての組合せで平均を
とり
2un ∑
{(xi − xn(1)(i+j) )T (xj − xn(2)(i+j) )}2
n(n − 1) i<j
n
Wn =
(4.2)
138
日本統計学会誌 第43巻 第1号 2013
を定義する．このとき，Wn は母集団分布に依らずに不偏性 E(Wn ) = tr(Σ2 ) をもつ．
Aoshima and Yata (2013a)，Yata and Aoshima (2013a) は一致性に次の定理を与えた．
定理 4.1
母集団分布に (A-iii) を仮定する．p → ∞，n → ∞ のとき，次が成り立つ．
(
Var
注意 4
Wn
tr(Σ2 )
)
=
4
{1 + o(1)} + O
n2
{
tr(Σ4 )
ntr(Σ2 )2
}
→ 0.
次の手順は，Wn の計算コストが O(pn2 ) のオーダーになり効率がよい．
(手順 1) (4.1) 式の xn(l)(k) , l = 1, 2 を各 k (= 3, . . . , 2n − 1) で計算する．
(手順 2) すべての i, j (1 ≤ i < j ≤ n) について手順 1 の xn(l)(i+j) を代入して un {(xi −
xn(1)(i+j) )T (xj − xn(2)(i+j) )}2 を計算し，それらの平均をとって Wn を得る．
定理 4.1 から分かるように，拡張クロスデータ行列法はクロスデータ行列法よりも漸近分
散が小さな不偏推定量を構築できる．Chen et al. (2010) は，U-統計量に基づいて tr(Σ2 )
∑n
∑n
の不偏推定量 {n(n − 1)}−1 i6=j (xTi xj )2 − 2{n(n − 1)(n − 2)}−1 i6=j6=k xTi xj xTj xk +
∑n
{n(n − 1)(n − 2)(n − 3)}−1 i6=j6=k6=l xTi xj xTk xl を与えた．これは Wn と同等な漸近分散
をもつが，計算コストが O(pn4 ) と非常に大きく，実用には向かない．拡張クロスデータ
行列法は，計算コストが O(pn2 ) であり，さらに，母集団に正規分布を仮定した場合には
Var{Wn /tr(Σ2 )} = 4{1 + o(1)}/n2 + 8tr(Σ4 ){1 + o(1)}/{ntr(Σ2 )2 } となり，これは正規
c2 ) と同等な漸近分散をもつ．
分布限定で提案された前掲の tr(Σ
4.2
高次元における相関ベクトルの検定
母集団に p + K 次元の分布を考え，n 個のデータ x1(∗) , . . . , xn(∗) を無作為に抽出した
とする．ただし，xj(∗) = (xTj , y1j , . . . , yKj )T ，j = 1, . . . , n とする．ここで，yij の分散
を Var(yij ) = σi2 ，xj と yij の相関ベクトルを Corr(xj , yij ) = ρi ，共分散ベクトルを
Cov(xj , yij ) = σ i とおく．そのとき，次の多重検定を考える．
Hi0 : ρi = 0 vs.
Hi1 : ρi 6= 0, i = 1, . . . , K.
(4.3)
高次元における相関の検定は，グラフィカルモデリングやパス解析等に用いられる．詳細
は，Wille et al. (2004)，Drton and Perlman (2007)，Hero and Rajaratnam (2011) 等や
Aoshima and Yata (2011a) を参照のこと．本節では，拡張クロスデータ行列法を使って検
定統計量を構築し，それに基づく多重検定が高い精度をもつことを示す．
各 i (= 1, . . . , K) で，k = 3, . . . , 2n − 1 について
y in(1)(k) = n−1
(1)
∑
j∈V
n(1)(k)
yij ,
y in(2)(k) = n−1
(2)
∑
j∈V
n(2)(k)
yij
139
高次元データの統計的方法論
と (4.1) の xn(1)(k) , xn(2)(k) を計算し，拡張クロスデータ行列法による検定統計量を
Tbn,σ i =
n
)(
)
)T (
)(
2un ∑ (
xt − xn(2)(s+t) yis − y in(1)(s+t) yit − y in(2)(s+t)
xs − xn(1)(s+t)
n(n − 1) s<t
(4.4)
で定義する．このとき，母集団分布に依らず不偏性 E(Tbn,σ i ) =k σ i k2 をもつ．次の 2 つを
仮定する．
(A-iv) p → ∞ のとき
tr(Σ4 )
→ 0.
tr(Σ2 )2
√
(A-v) 各 i (= 1, . . . , K) で，p → ∞, n → ∞, k σ i k6= 0 のとき lim inf
σi2
tr(Σ2 )
n k σ i k2
> 0.
(A-iv) は，“p → ∞ のとき λ1 /tr(Σ2 )1/2 → 0”という，固有値に関する条件と同値である．
このとき，Yata and Aoshima (2013a) は次の定理を与えた．
定理 4.2
母集団分布に (A-ii) を仮定する．(A-iv)–(A-v) と，yij に関するある正則条件
のもとで，p → ∞，n → ∞ のとき次が成り立つ．
Tbn,σ i − k σ i k2
√
⇒ N (0, 1),
Sin(y) 2Wn /n
i = 1, . . . , K.
ここで，Sin(y) は yi1 , . . . , yin の不偏標本分散，Wn は (4.2) 式で与えられるものである．
D = {i | ρi 6= 0, i = 1, . . . , K}, D c = {i | ρi = 0, i = 1, . . . , K} とおき，zα/K を
N (0, 1) の上側 α/K 点とする．(4.3) に関する検定方式を α ∈ (0, 1/2) に対して
Hi0 を棄却 ⇐⇒
Tbn,σ i
√
> zα/K ,
Sin(y) 2Wn /n
i = 1, . . . , K
(4.5)
b = {i | Hi0 を棄却, i = 1, . . . , K} とする．定理 4.2 から，family-wise error
で定義し，D
b 6= ∅) と検出力 P (D ⊂ D)
b について次の系が成り立つ．
rate (FWER) である P (D c ∩ D
系 4.1
母集団分布に (A-ii) を仮定する．(A-iv)–(A-v) と，yij に関するある正則条件の
もとで，p → ∞，n → ∞ のとき検定方式 (4.5) について次が成り立つ．
(
)
2
∑
n
k
σ
k
i
c
b 6= ∅) ≤ α, lim inf P (D ⊂ D)
b ≥1−
Φ zα/K − √
lim sup P (D ∩ D
.
σi2 2tr(Σ2 )
i∈D
注意 5
√
定理 4.2 から，漸近 P 値である Pi = 1−Φ{Tbn,σ i /(Sin(y) 2Wn /n)}, i = 1, . . . , K
を用いれば，Holm 法 (Holm (1979)) や false discovery rate (FDR) を調整する B-H 法
(Benjamini and Hochberg (1995)) も考えられる．
140
日本統計学会誌 第43巻 第1号 2013
FWER は保守的で，FDR 等に比べ検出力が小さい傾向にある．そこで，Yata and Aoshima
(2013a) は，事前に FWER と検出力の目標精度を設定し，これらを漸近的に満足するよう
な検定方式を構築した．与えられる α, β ∈ (0, 1/2), ∆L = o{tr(Σ2 )1/2 /tr(Σ)} (> 0) に対
b 6= ∅) ≤ α を満たし，min k σ i k2 /{σ 2 tr(Σ)} ≥ ∆L のときの検
して，FWER が P (D c ∩ D
i
i∈D
b ≥ 1 − β となるような検定方式を考える．ただし，k σ i k2 /{σ 2 tr(Σ)}
出力が P (D ⊂ D)
i
は yij から xj への寄与率を表す (注意 7 参照)．系 4.1 で与えた検出力から，標本数 n が
n≥
(zα/K + zβ/K )
∆L tr(Σ)
√
2tr(Σ2 ) (= C とおく)
b ≥ 1 − β を満たす．Yata and
を満たすならば，検定方式 (4.5) は漸近的に P (D ⊂ D)
Aoshima (2013a) は，未知パラメータ Σ に依存する C を初期データで推定し，要求される
精度を漸近的に保証する 2 段階検定法を与えた．
[2 段階検定法]
(手順 1) 初期標本数 m (≥ 4) を lim inf p→∞ m/C > 0 と lim supp→∞ m/C < ∞ を満たす
ように選ぶ．無作為に m 個のデータ xj(∗) , j = 1, . . . , m を抽出し，xj , j = 1, . . . , m
の標本共分散行列 S m と (4.2) 式に従って Wm を計算する．そのとき，全標本数を
m}
{
l (z
√
α/K + zβ/K )
2Wm + 2zβ/K (zα/K + zβ/K )
N = max m,
∆L tr(S m )
で定義する．
(手順 2) 無作為に N − m 個の追加データ xj(∗) , j = m + 1, . . . , N を抽出し，xj , j =
1, . . . , N の標本共分散行列 S N と各 i (= 1, . . . , K) の SiN (y) を計算し，(4.4) 式に
従って TbN,σ i を計算する．そのとき，各 i (= 1, . . . , K) で検定方式を次のように定
義する．
{1 − 2zβ/K (zα/K + zβ/K )/N }SiN (y) tr(S N )∆L zα/K
.
Hi0 を棄却 ⇐⇒ TbN,σ i >
zα/K + zβ/K
(4.6)
Yata and Aoshima (2013a) は次の定理を与えた．
定理 4.3
母集団分布に (A-ii) を仮定する．(A-iv) と，yij に関するある正則条件のもと
で，p → ∞ のとき検定方式 (4.6) について次が成り立つ．
b 6= ∅) ≤ α,
lim sup P (D c ∩ D
2
b ≥ 1 − β when min k σ i k ≥ ∆L .
lim inf P (D ⊂ D)
2
i∈D σi tr(Σ)
高次元データの統計的方法論
注意 6
141
lim inf p→∞ p∆L > 0 ならば，定理 4.3 は N/p = op (1) なる高次元小標本
の枠組みでも成り立つことに注意する．Aoshima and Yata (2010) は，各種 2 段階推
測法における初期標本数 m の決め方と推測の 2 次漸近一致性を研究した．上記の 2 段
階検定法の場合，tr(Σ2 )1/2 /tr(Σ) ≥ p−1/2 が成立することに注意して，初期標本数を
√
√
m = max{4, d(zα/K + zβ/K ) 2/(∆L p) + 2zβ/K (zα/K + zβ/K )e} で決める．そのと
き，tr(Σ2 )1/2 /tr(Σ) = O(p−1/2 ) なる仮定のもとで定理 4.3 が成り立つ．
注意 7
k σ i k2 /{σi2 tr(Σ)} ∈ [0, 1] は回帰モデルの重相関係数になっており，yij から xj
への寄与率を表す．詳細は，Yata and Aoshima (2013a) の Remark 8 を参照されたい．事
前に与える ∆L は，各 yij から xj への寄与率に関して，有意と考えられる下限の設定値で
ある．2 段階検定法は，∆L 以上の寄与率を有する変数 yij (i ∈ D) をすべて選択する確率
が，漸近的に 1 − β 以上になるような検定方式である．
4.3
実データ解析
本節では，2 段階検定法の実例として，Yata and Aoshima (2013a) で扱ったマイクロアレ
イデータのパス解析の一部を紹介する．Wille et al. (2004) のシロイヌナズナのマイクロア
レイデータを用いた．これは，39 (= K) 個のイソプレノイドに関する遺伝子と 795 (= p) 個
の追加遺伝子からなる 834 次元データである．推測の精度を α = 0.05，β = 0.1，∆L = 0.1
b 6= ∅) ≤ 0.05，且つ，39 個の遺伝子から寄与率
と設定した．すなわち，FWER P (D c ∩ D
が 10% 以上の遺伝子をすべて選択する確率が 90% 以上になるように検定を設計した．
まず，注意 6 から初期標本数 m を
√
m}
{ l (z
α/39 + zβ/39 ) 2
m = max 4,
+ 2zβ/39 (zα/39 + zβ/39 ) = 36
√
0.1 p
とした．36 個の初期データを抽出し，tr(S m ) = 440 と Wm = 13935 を得た．そのとき，
全標本数は
{
N = max m,
l (z
√
m}
+ zβ/39 ) 2 √
13935 + 2zβ/39 (zα/39 + zβ/39 ) = 55
0.1 × 440
α/39
となった．N − m = 19 個の追加データを抽出して S N , SiN (y) , TbN,σ i , i = 1, . . . , 39 を計
算し，検定方式 (4.6) により 39 個の遺伝子を検定した．その結果，設定した推測の精度で，
b = {1, . . . , 39} \ {6, 7, 13, 14, 15, 16, 17, 20} に該当する 31 個の遺伝子は寄与率が高いと
D
b のパス解析を行い，これらが予
判定された．Yata and Aoshima (2013a) は，遺伝子群 D
測の意味で有効に選ばれていることを確認している．
142
5.
日本統計学会誌 第43巻 第1号 2013
高次元の平均ベクトルに関する推定と検定
本節では，母集団が k 個あると想定し，各母集団 (πi ) は平均に p 次のベクトル µi ，共
分散行列に p 次の非負定値対称行列 Σi (≥ O) をもつと仮定する．母集団 πi , i = 1, . . . , k
の高次元データに次のモデルを考える．
xij = Γi wij + µi ,
j = 1, . . . , ni (≥ 4).
ここで，Γi は Γi ΓTi = Σi なる p × ri 行列，wij は E(w ij ) = 0, Var(wij ) = I ri と
し，wij = (wi1j , . . . , wiri j )T の各成分は 4 次モーメントが一様有界と仮定する．母集団
πi , i = 1, . . . , k の分布には，次を仮定する．
(A-vi)
2
2
E(wiqj
wisj
) = 1，E(wiqj wisj witj wiuj ) = 0，q 6= s, t, u.
Σi の最大固有値を λi1 とおく．各母集団の共分散行列に次を仮定する．
λi1
→ 0, p → ∞; i = 1, . . . , k.
tr(Σ2i )1/2
∑k
平均ベクトルの 1 次結合 µ = i=1 bi µi に関する推測を考える．次元数 p が固定の枠組
(A-vii)
みでは，Aoshima et al. (2002)，Aoshima and Takada (2004)，Aoshima and Yata (2010)
等によって，事前に設定された推測の精度を満たすための各種 2 段階推測法が研究された．
次元数が p → ∞ のとき，k = 2 且つ (b1 , b2 ) = (1, −1) なる高次元データの 2 標本問題につ
いては，Bai and Saranadasa (1996)，Srivastava (2007)，Chen and Qin (2010)，Aoshima
and Yata (2011a, 2013b) の研究がある．特に，Aoshima and Yata (2011a, 2013b) は，母
集団分布と共分散行列に関する仮定を上記の (A-vi)–(A-vii) にまで緩め，より洗練された
解を与えている．
∑k
各母集団 πi から無作為に ni 個のデータを抽出して，T n =
i=1 bi xini を定義する．
∑ni
ただし，n = (n1 , . . . , nk )，xini =
j=1 xij /ni である．各母集団 πi の標本共分散行列
∑
n
i
b n = ∑k b2 tr(S in )/ni とお
S ini = (ni −1)−1 j=1
(xij −xini )(xij −xini )T に基づいて，Σ
i
i=1 i
∑
k
2
2
2
4
b n ) = 0，Var(k T n −µ k −Σ
b n) = 2
く．このとき，E(k T n −µ k −Σ
i=1 bi tr(Σi )/{ni (ni −
∑k
1)} + 4 i<j b2i b2j tr(Σi Σj )/(ni nj ) (= δ とおく) である．T n のノルムについて，Aoshima
and Yata (2011a, 2013b) は次の定理を与えた．
定理 5.1
(A-vi)–(A-vii) のもとで，p → ∞, ni → ∞, i = 1, . . . , k のとき次が成り立つ．
bn
k T n − µ k2 −Σ
⇒ N (0, 1).
1/2
δ
b n とおく．このとき，E(Tbn,µ ) =k µ k2 である．Aoshima and Yata
Tbn,µ =k T n k2 −Σ
(2011a, 2013b) は次の定理を与えた．
高次元データの統計的方法論
δ −1
定理 5.2
∑k
143
µT Σi µ/ni = o(1) を仮定する．(A-vi)–(A-vii) のもとで，p → ∞,
i=1
ni → ∞, i = 1, . . . , k のとき次が成り立つ．
Tbn,µ − k µ k2
⇒ N (0, 1).
δ 1/2
未知の δ は次のように推定する．
δ̂ = 2
k
∑
b4i Wini /{ni (ni − 1)} + 4
i=1
k
∑
b2i b2j tr(S ini S jnj )/(ni nj ).
i<j
ただし，Wini は (4.2) 式に従って計算する tr(Σ2i ) の不偏推定量である．次の系を得る．
系 5.1
(A-vi)–(A-vii) のもとで，p → ∞, ni → ∞, i = 1, . . . , k のとき次が成り立つ．
bn
k T n − µ k2 −Σ
さらに，δ −1
∑k
i=1
δ̂ 1/2
⇒ N (0, 1).
µT Σi µ/ni = o(1) を仮定すれば次が成り立つ．
Tbn,µ − k µ k2
δ̂ 1/2
⇒ N (0, 1).
これらの結果から，平均ベクトルに関する推測に幾つかの応用が考えられる．Aoshima
and Yata (2013b) は，
b n − zα/2 δ̂ 1/2 , 0} ≤k T n − µ k 2 ≤ Σ
b n + zα/2 δ̂ 1/2 }
Rµ = {µ ∈ Rp : max{Σ
なる µ の信頼領域を考えた．ここで，zα は N (0, 1) の上側 α 点を表す．系 5.1 から，(A-
vi)–(A-vii) のもとで，p → ∞, ni → ∞, i = 1, . . . , k のとき次が成り立つ．
P (µ ∈ Rµ ) = 1 − α + o(1).
また，Aoshima and Yata (2013b) は，k = 2 且つ (b1 , b2 ) = (1, −1) のとき 2 標本問題とな
るような
H0 : µ = 0 vs. H1 : µ 6= 0
なる µ の検定を考え，α ∈ (0, 1/2) に対して検定方式を
H0 を棄却 ⇐⇒
で定義した．系 5.1 から，δ −1
∑2
i=1
Tbn,µ
δ̂ 1/2
> zα
µT Σi µ/ni = o(1) を仮定すれば，(A-vi)–(A-vii) のも
とで，p → ∞，ni → ∞，i = 1, 2 のとき次が成り立つ．
)
(
k µ k2
−
z
Size = α + o(1), Power = Φ
α + o(1).
δ 1/2
Aoshima and Yata (2013b)，Yata and Aoshima (2012b) は，上記以外にも，高次元の平
均ベクトルに関する各種推測を扱っている．事前に設定される検出力の目標値を達成する
ための 2 段階推測法については，Aoshima and Yata (2011a, b) を参照のこと．
144
日本統計学会誌 第43巻 第1号 2013
高次元データの判別分析
6.
高次元データの 2 群の判別分析を考える．標本共分散行列 S ini の逆行列が存在しないた
め，Fisher の線形判別方式や 2 次判別方式は使えない．2 群の共分散行列が等しいと仮定
すれば，Dudoit et al. (2002) や Bickel and Levina (2004) による標本共分散行列の対角成
分だけを使った判別方式，Srivastava and Kubokawa (2007) のリッジ型逆行列による判別
方式，Yata and Aoshima (2012a) のノイズ掃き出し法を用いた共分散行列の逆行列推定
による判別方式がある．しかし，2 群の共分散行列が等しいと仮定する問題設定の単純化
は，高次元データが本来もつ 2 群の差異に関する情報に目を瞑ることになる．2 群の共分散
行列が等しいことを仮定しない場合，Dudoit et al. (2002) による標本共分散行列の対角成
分だけを使った判別方式，Hall et al. (2008)，Chan and Hall (2009)，Aoshima and Yata
(2013a) 等のユークリッド距離に基づく判別方式，Aoshima and Yata (2011a) による高次
元小標本の幾何学的表現に基づく判別方式がある．他にも，Aoshima and Yata (2011a)，
Fan and Fan (2008) 等による変数選択に基づく判別方式がある．
2 つの母集団 πi , i = 1, 2 は，平均に p 次のベクトル µi ，共分散行列に p 次の正定値対
称行列 Σi (> O) をもつとする．各母集団 πi から，ni (≥ 4) 個のデータ xi1 , . . . , xini を
無作為に抽出する．判別対象のデータを x0 とする．判別の精度について，π1 の x0 を π2
に誤判別する確率を e(2|1)，π2 の x0 を π1 に誤判別する確率を e(1|2) で表す．本節では，
Aoshima and Yata (2013a) によるユークリッド距離に基づく判別方式と，Aoshima and
Yata (2011a) による高次元小標本の幾何学的表現に基づく判別方式を紹介し，他の判別方
式も含め性能を数値的に比較する．
6.1
ユークリッド距離に基づく判別方式
Aoshima and Yata (2013a) は，ユークリッド距離に基づく次のような判別関数を考えた．
(
x1n1 + x2n2 )T
tr(S 1n1 ) tr(S 2n2 )
ωa (x0 ) = x0 −
(x2n2 − x1n1 ) −
+
.
(6.1)
2
2n1
2n2
判別方式は，ωa (x0 ) < 0 のとき x0 ∈ π1 ，ωa (x0 ) ≥ 0 のとき x0 ∈ π2 とする．∆ =k µ1 −µ2 k2
とおき，次の 2 つを仮定する．
(A-viii)
(A-ix)
(µ1 − µ2 )T Σi (µ1 − µ2 )
→ 0, p → ∞, i = 1, 2.
∆2
maxj=1,2 tr(Σ2j )
→ 0, p → ∞, i = 1, 2.
ni ∆ 2
Aoshima and Yata (2013a) は次の定理を与えた．
定理 6.1
(A-viii)–(A-ix) のもとで，p → ∞ のとき判別方式 (6.1) について次が成り
立つ．
e(2|1) → 0,
e(1|2) → 0.
145
高次元データの統計的方法論
上記の結果は，母集団 πi , i = 1, 2 の分布に依存しない．つまり，判別方式 (6.1) は，母集団
分布に (A-vi) 等を仮定せずに，平均間のユークリッド距離 ∆ を使って精度よく判別している．
もしも，母集団分布に (A-vi) が仮定できるならば，判別関数 ωa (x0 ) の漸近正規性が次のように
∑2
得られる．いま，κi = tr(Σ2i )/ni +tr(Σ1 Σ2 )/ni0 + j=1 tr(Σ2j )/{2nj (nj −1)}, i(6= i0 ) = 1, 2
とおく．次を仮定する．
(A-x)
(µ1 − µ2 )T Σi (µ1 − µ2 )
= o(1), i = 1, 2.
κi
ここで，仮定 (A-x) のもと Var(ωa (x0 )|x0 ∈ πi )/κi → 1, i = 1, 2 である．Aoshima and
Yata (2013a) は次の定理を与えた．
定理 6.2
母集団分布に (A-vi) を仮定する．(A-vii) と (A-x) と適当な正則条件のもと
で，p → ∞, ni → ∞, i = 1, 2 のとき，判別関数 ωa (x0 ) について次が成り立つ．
ωa (x0 ) − (−1)i ∆/2
1/2
κi
⇒ N (0, 1),
x0 ∈ πi , i = 1, 2.
Aoshima and Yata (2013a) は，誤判別確率に対して事前に設定される上限を超えないよ
うな精度保証を考え，判別方式を上記の漸近正規性に基づいて構築した．さらに，母集団
の個数が 3 個以上の場合にも，同様の精度を保証する判別方式を開発した．
6.2
幾何学的表現に基づく 2 次判別方式
Aoshima and Yata (2011a) は，高次元小標本の幾何学的表現に基づいて次のような判別
関数を考えた．
ωb (x0 ) =
p k x0 − x1n1 k2
p k x0 − x2n2 k 2
−
− p log
tr(S 1n1 )
tr(S 2n2 )
{
tr(S 2n2 )
tr(S 1n1 )
}
−
p
p
+ .
n1
n2
(6.2)
判別方式は，ωb (x0 ) < 0 のとき x0 ∈ π1 ，ωb (x0 ) ≥ 0 のとき x0 ∈ π2 とする．ここで，
母集団分布に (A-vi) を仮定し，Σi , i = 1, 2 の固有値に球形条件 (3.1) を仮定すれば，
x0 ∈ πi (i = 1, 2) について次が成り立つ．
k x0 − µi k2
= 1 + op (1),
tr(Σi )
p → ∞.
これは，規準化したデータの分布の半径が一定値に近づくことを意味する．いま，γi =
{tr(Σ1 ) − tr(Σ2 )}/tr(Σi ), i = 1, 2 とおき，limp→∞ γ2 → c (6= 0) となる場合を考える．幾
つかの正則条件のもとで，p → ∞ のとき x0 ∈ π1 について次が成り立つ．
ωb (x0 )
∆
= −γ2 + log (1 + γ2 ) −
+ oP (1) ≤ −c + log(1 + c) + oP (1) < 0.
p
tr(Σ2 )
同様に，x0 ∈ π2 に対しては漸近的に正となる．すなわち，limp→∞ γ2 → c (6= 0) なる意味
で Σ1 と Σ2 に差異がある場合は，判別方式 (6.2) の誤判別確率は漸近的に 0 となる．
146
日本統計学会誌 第43巻 第1号 2013
次に，Σ1 と Σ2 に十分な差異がなく limp→∞ γ2 → 0 となる場合を考える．その場合，
p → ∞ のときの ωb (x0 )/p の主要項を
γi − log (1 + γi ) +
∆ + {tr(Σ1 ) − tr(Σ2 )}2 /{2tr(Σi )}
∆
=
{1 + o(1)},
tr(Σi )
tr(Σi )
i = 1, 2
と表せるので，∆? = ∆ + mini=1,2 {tr(Σ1 ) − tr(Σ2 )}2 /{2tr(Σi )} は 2 群間の差異を表すパ
ラメータとみなせる．次の 2 つを仮定する．
(A-xi)
(A-xii)
(µ1 − µ2 )T Σi (µ1 − µ2 )
tr(Σ2i ){tr(Σ1 ) − tr(Σ2 )}2
→
0,
→ 0, p → ∞, i = 1, 2.
∆2?
tr(Σi )2 ∆2?
maxj=1,2 tr(Σ2j )
→ 0, p → ∞, i = 1, 2.
ni ∆2?
Aoshima and Yata (2011a, 2013a) は次の定理を与えた．
定理 6.3
母集団分布に (A-vi) を仮定する．(A-xi) と (A-xii) と適当な正則条件のもと
で，p → ∞ のとき判別方式 (6.2) について次が成り立つ．
e(2|1) → 0,
e(1|2) → 0.
Aoshima and Yata (2011a) は，判別関数 ωb (x0 ) の高次元漸近正規性を示し，誤判別確
率に対して事前に設定される上限を超えないような判別方式を与えた．
6.3
判別方式の性能比較
判別方式 (6.2) は，∆? を通して，µi , i = 1, 2 の差異 ∆ だけでなく Σi , i = 1, 2 の差異も
考慮することで，判別精度を向上させている．∆? が大きいほど誤判別確率のゼロへの収束
は速い．ただし，母集団分布に (A-vi) の仮定を必要とする．一方で，判別方式 (6.1) は，∆
の大きさだけで判別を行う．母集団分布に仮定を必要とせず，簡易で頑健な方法といえる．
判別方式 (6.1) と (6.2) について，他の判別方式の検討も含め，数値的に性能を比較し
た．次のようにシミュレーション実験を設計した．2 つの母集団の平均ベクトルは µ1 =
0, µ2 = (p−1/6 , . . . , p−1/6 )T とした．2 つの母集団の共分散行列は Σ1 = 0.8B(0.3|i−j| )B ，
1/3
Σ2 = 1.2B(0.3|i−j| )B ，B = diag[{0.5 + 1/(p + 1)}1/2 , . . . , {0.5 + p/(p + 1)}1/2 ] とした．
1/3
そのとき，∆ =k µ1 − µ2 k2 = p2/3 ，∆? = p/15 + p2/3 となる．任意の自然数 ni , i = 1, 2
について，(A-viii) から (A-xii) が成り立つ．標本数は n1 = 10, n2 = 20 とした．母集団
πi , i = 1, 2 の分布に 2 つの場合を考えた．
(a) p 次元正規分布 Np (µi , Σi ), p = 2s , s = 5, . . . , 11.
(b) 平均ベクトルが µi ，共分散行列が Σi ，自由度 ν = 10(10)60 の p = 500 次元 t 分布．
高次元データの統計的方法論
(a) p 次元の正規分布, p = 2s , s = 5, . . . , 11.
図7
147
(b) 500 次元の自由度 ν = 10(10)60 の t 分布.
A: ωa (x0 )，B: ωb (x0 )，C: DLDA，D: DQDA，E: HM-LSVM による誤判別確率．
これら 2 つの場合は，(b) の t 分布の自由度 ν が大きくなると，(a) の正規分布 Np (µi , Σi )
に近づくという関係にある．
判別方式として (6.1) と (6.2) の他に，Dudoit et al. (2002) による標本共分散行列の対角
成分だけを使った線形判別方式 (DLDA)，2 次判別方式 (DQDA)，そして，Vapnic (1999)
によるハードマージン・線形サポートベクターマシン (HM-LSVM) も考え，合計 5 つの判
別方式を比較した．サポートベクターマシンに HM-LSVM を用いた理由は，p > n1 + n2
なる高次元小標本の設定により，適当な超平面で完全分離可能だからである．判別対象が
x0 ∈ π1 と x0 ∈ π2 のそれぞれの場合で，上記の 5 つの判別方式が x0 を正しく判別するか
を確認した．これを独立に 2000 回繰り返し，x0 ∈ π1 と x0 ∈ π2 のそれぞれの場合におけ
る誤判別の割合 e(2|1), e(1|2) を計算し，平均誤判別確率として e = (e(2|1) + e(1|2))/2 を
図 7 にプロットした．このとき，e の標準偏差は 0.011 以下である．
図 7 の (a) から，(A-vi) を満たすときには，判別関数 ωb (x0 ) が最も性能が良いように見
受けられる．µi , i = 1, 2 だけでなく，Σi , i = 1, 2 の差異も考慮した効果が表れていると
いえよう．一方，図 7 の (b) を見ると，自由度 ν が小さく，そのため (A-vi) を満たさない
ときには，判別方式 (6.1) の性能が (6.2) の性能に勝ることが分かる．予想される通り，ν
が増加して (A-vi) を満たすようになれば，判別方式 (6.2) の性能が一番良くなる．DQDA
は高次元データに対して非常に大きなバイアスが生じるために，精度が著しく悪くなって
いる．判別方式 (6.1) と DLDA は常に安定した結果を与え，非常に頑健であるといえる．
特に，判別方式 (6.1) は，判別関数が定理 6.2 の漸近正規性を有し，それに基づいて誤判別
確率の調整ができる方法論であるため，非常に有用であるといえる．ここでは割愛するが，
シミュレーションの設定を変えても同様な結果を得ている．
謝辞
本研究は，科学研究費補助金 基盤研究 (B) 22300094 研究代表者: 青嶋 誠「高次元データ
148
日本統計学会誌 第43巻 第1号 2013
の理論と方法論の総合的研究」
，および，学術研究助成基金助成金 挑戦的萌芽研究 23650142
研究代表者: 青嶋 誠「高速で頑健かつ高精度な多変量統計手法の新展開」，若手研究 (B)
23740066 研究代表者: 矢田 和善「高次元小標本の理論的体系の構築」から研究助成を受け
ている．初稿を改訂するにあたって，編集委員の方から適切なコメントを頂いた．ここに
謝意を表したい．
参 考 文 献
Ahn, J., Marron, J. S., Muller, K. M. and Chi, Y.-Y. (2007). The high-dimension, low-sample-size geometric
representation holds under mild conditions, Biometrika, 94, 760–766.
Alon, U., Barkai, N., Notterman, D. A., Gish, K., Ybarra, S., Mack, D. and Levine, A. J. (1999). Broad patterns
of gene expression revealed by clustering analysis of tumor and normal colon tissues probed by oligonucleotide
arrays, Proc. Natl. Acad. Sci. USA, 96, 6745–6750.
Aoshima, M. and Takada, Y. (2004). Asymptotic second-order efficiency for multivariate two-stage estimation
of a linear function of normal mean vectors, Seq. Anal., 23, 333–353.
Aoshima, M. and Yata, K. (2010). Asymptotic second-order consistency for two-stage estimation methodologies
and its applications, Ann. Inst. Statist. Math., 62, 571–600.
Aoshima, M. and Yata, K. (2011a). Two-stage procedures for high-dimensional data, Seq. Anal. (Editor’s special
invited paper ), 30, 356–399.
Aoshima, M. and Yata, K. (2011b). Authors’ response, Seq. Anal., 30, 432–440.
Aoshima, M. and Yata, K. (2011c). Effective methodologies for statistical inference on microarray studies, in
Prostate Cancer—From Bench to Bedside, edited by Spiess, P. E., pp. 13–32, InTech.
Aoshima, M. and Yata, K. (2013a). A distance-based, misclassification rate adjusted classifier for multiclass,
high-dimensional data, Ann. Inst. Statist. Math., revised.
Aoshima, M. and Yata, K. (2013b). Asymptotic normality for inference on multisample, high-dimensional mean
vectors under mild conditions, Methodol. Comput. Appl. Probab., doi:10.1007/s11009-013-9370-7, in press.
青嶋誠・矢田和善 (2013). 「論説：高次元小標本における統計的推測」『数学』65, 225–247．
Aoshima, M., Takada, Y. and Srivastava, M. S. (2002). A two-stage procedure for estimating a linear function
of k multinormal mean vectors when covariance matrices are unknown, J. Statist. Plann. Inference, 100,
109–119.
Bai, Z. and Saranadasa, H. (1996). Effect of high dimension: by an example of a two sample problem, Statist.
Sinica, 6, 311–329.
Baik, J. and Silverstein, J. W. (2006). Eigenvalues of large sample covariance matrices of spiked population
models, J. Multivar. Anal., 97, 1382–1408.
Baik, J., Ben Arous, G. and Péché, S. (2005). Phase transition of the largest eigenvalue for nonnull complex
sample covariance matrices, Ann. Probab., 33, 1643–1697.
Benjamini, Y. and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach
to multiple testing, J. R. Stat. Soc. Ser. B Stat. Methodol., 57, 289–300.
Bickel, P. J. and Levina, E. (2004). Some theory for Fisher’s linear discriminant function, ‘naive Bayes’, and
some alternatives when there are many more variables than observations, Bernoulli , 10, 989–1010.
Chan, Y.-B. and Hall, P. (2009). Scale adjustments for classifiers in high-dimensional, low sample size settings,
Biometrika, 96, 469–478.
Chen, S. X. and Qin, Y.-L. (2010). A two-sample test for high-dimensional data with applications to gene-set
testing, Ann. Statist., 38, 808–835.
Chen, S. X., Zhang, L.-X. and Zhong, P.-S. (2010). Tests for high-dimensional covariance matrices, J. Am.
Statist. Assoc., 105, 810–819.
Drton, M. and Perlman, M. D. (2007). Multiple testing and error control in Gaussian graphical model selection,
Statist. Sci., 22, 430–449.
高次元データの統計的方法論
149
Dudoit, S., Fridlyand, J. and Speed, T. P. (2002). Comparison of discrimination methods for the classification
of tumors using gene expression data, J. Am. Statist. Assoc., 97, 77–87.
Fan, J. and Fan, Y. (2008). High-dimensional classification using features annealed independence rules, Ann.
Statist., 36, 2605–2637.
Golub, T. R., Slonim, D. K., Tamayo, P., Huard, C., Gaasenbeek, M., Mesirov, J. P., Coller, H., Loh, M. L.,
Downing, J. R., Caligiuri, M. A., Bloomfield, C. D. and Lander, E. S. (1999). Molecular classification of
cancer: class discovery and class prediction by gene expression monitoring, Science, 286, 531–537.
Hall, P., Marron, J. S. and Neeman, A. (2005). Geometric representation of high dimension, low sample size
data, J. R. Stat. Soc. Ser. B Stat. Methodol., 67, 427–444.
Hall, P., Pittelkow, Y. and Ghosh, M. (2008). Theoretical measures of relative performance of classifiers for
high dimensional data with small sample sizes, J. R. Stat. Soc. Ser. B Stat. Methodol., 70, 159–173.
Hero, A. and Rajaratnam, B. (2011). Large-scale correlation screening, J. Am. Statist. Assoc., 106, 1540–1552.
Holm, S. (1979). A simple sequentially rejective multiple test procedure, Scand. J. Statist., 6, 65–70.
Johnstone, I. M. (2001). On the distribution of the largest eigenvalue in principal components analysis, Ann.
Statist., 29, 295–327.
Johnstone, I. M. and Lu, A. Y. (2009). On consistency and sparsity for principal components analysis in high
dimensions, J. Am. Statist. Assoc., 104, 682–693.
Jung, S. and Marron, J. S. (2009). PCA consistency in high dimension, low sample size context, Ann. Statist.,
37, 4104–4130.
Jung, S., Sen, A. and Marron, J. S. (2012). Boundary behavior in high dimension, low sample size asymptotics
of PCA, J. Multivar. Anal., 109, 190–203.
Lee, S., Zou, F. and Wright, F. A. (2010). Convergence and prediction of principal component scores in highdimensional settings, Ann. Statist., 38, 3605–3629.
Paul, D. (2007). Asymptotics of sample eigenstructure for a large dimensional spiked covariance model, Statist.
Sinica, 17, 1617-1642.
Singh, D., Febbo, P. G., Ross, K., Jackson, D. G., Manola, J., Ladd, C., Tamayo, P., Renshaw, A. A., D’Amico,
A. V., Richie, J. P., Lander, E. S., Loda, M., Kantoff, P. W., Golub, T. R. and Sellers, W. R. (2002). Gene
expression correlates of clinical prostate cancer behavior, Cancer Cell , 1, 203–209.
Srivastava, M. S. (2005). Some tests concerning the covariance matrix in high dimensional data, J. Jpn. Statist.
Soc., 35, 251–272.
Srivastava, M. S. (2007). Multivariate theory for analyzing high dimensional data, J. Jpn. Statist. Soc., 37,
53–86.
Srivastava, M. S. and Kubokawa, T. (2007). Comparison of discrimination methods for high dimensionnal data,
J. Jpn. Statist. Soc., 37, 123–134.
Vapnic, V. N. (1999). The Nature of Statistical Learning Theory (second ed.), Springer-Verlag, New York.
Wille, A., Zimmermann, P., Vranová, E., Fürholz, A., Laule, O., Bleuler, S., Hennig, L., Prelic, A., von Rohr,
P., Thiele, L., Zitzler, E., Gruissem, W. and Bühlmann, P. (2004). Sparse graphical Gaussian modeling of
the isoprenoid gene network in Arabidopsis thaliana, Genome Biol., 5, R92.
Yata, K. (2010). Effective two-stage estimation for a linear function of high-dimensional Gaussian means, Seq.
Anal., 29, 463–482.
Yata, K. and Aoshima, M. (2009). PCA consistency for non-Gaussian data in high dimension, low sample size
context, Comm. Statist. Theory Methods, Special Issue Honoring Zacks, S. (ed. Mukhopadhyay, N.), 38,
2634–2652.
Yata, K. and Aoshima, M. (2010a). Effective PCA for high-dimension, low-sample-size data with singular value
decomposition of cross data matrix, J. Multivar. Anal., 101, 2060–2077.
Yata, K. and Aoshima, M. (2010b). Intrinsic dimensionality estimation of high-dimension, low sample size
data with d-asymptotics, in Comm. Statist. Theory Methods, Special Issue Honoring M. Akahira, edited by
Aoshima, M., 39, 1511–1521.
Yata, K. and Aoshima, M. (2012a). Effective PCA for high-dimension, low-sample-size data with noise reduction
via geometric representations, J. Multivar. Anal., 105, 193–215.
150
日本統計学会誌 第43巻 第1号 2013
Yata, K. and Aoshima, M. (2012b). Inference on high-dimensional mean vectors with fewer observations than
the dimension, Methodol. Comput. Appl. Probab., 14, 459–476.
Yata, K. and Aoshima, M. (2013a). Correlation tests for high-dimensional data using extended cross-datamatrix methodology, J. Multivar. Anal., 117, 313–331.
Yata, K. and Aoshima, M. (2013b). PCA consistency for the power spiked model in high-dimensional settings,
J. Multivar. Anal., doi:10.1016/j.jmva.2013.08.003, in press.