高次元現象の統計数理
青嶋誠(筑波大数理物質)∗

1. はじめにゲノム科学・情報工学・金融工学などの現代科学の一つの特徴は，データがもつ次元数の膨大さにある．例えば，DNA マイクロアレイによる遺伝子発現データを見ると，次元数が数万にものぼる一方で，標本数は 100 にも満たないという事例が多い．高次元
データは，高次元になるほど標本を入手することが難しく，一般にデータの次元数 p と標本数 n にはp ≫ n (もしくは，p > n)といった大小関係がある．これを強調して，高次元小標本データとよぶこともある．
2005 年頃から，次世代シーケンサーによる超高次元ゲノムデータも統計解析の対象となり，例えば Wang et al. [13] は，次元数が 300 万を超え，標本数は 40 というゲノムデータを扱っている．近年，高次元小標本データの解析は，様々な分野で需要が高まっている．本講演は，次世代シーケンサーによる超高次元ゲノムデータをはじめ，天文
宇宙学における天体の分光マッピングによる波長データ ([11]) や，材料科学における
結晶構造データなど，高次元小標本データの解析例を幾つか紹介する．
高次元データは，豊富な情報を有するものの，それが巨大なノイズに埋もれている
ために見つけ難い．伝統的な多変量解析法では，高次元データの解析に精度を保証することができない．次元数が標本数を超えると，次元の呪いとして様々な高次元現象が現れ，多変量解析の理論と方法論の成立を阻むのである．特に非スパース性という現象に陥ると，データの潜在情報は高次元において膨張するノイズに埋没してしまい，
推測の強不一致という最悪な結果を招くことにもなる．
高次元データの解析には，高次元データ特有の理論と方法論が必要になる．著者の研究グループは，高次元統計解析と非スパースモデリングという新しい理論と方法論を考案した．本講演は，様々な高次元現象を説明し，たった数十程度の小標本であっても高速かつ高精度に解析するための，高次元現象の統計数理を紹介する．

2. 高次元データの非スパース性
1990 年代後半，高次元データを解析するための方法論として，機械学習が目覚ましく発展した．スタンフォード大学の Donoho 教授や Tibshirani 教授らは，スパースモデリ
ングを発展させた．これは，高次元データの母共分散行列にスパース性を仮定して，高次元データの特徴を少数の変数や標本で捉えようとするもので，L1 正則化法や圧縮センシングがよく知られる ([7])．近年，多層ニューラルネットワークを用いた深層学習が盛んに研究されている．これもスパース性を仮定して，データは大標本であること
本研究は，筑波大学の矢田和善准教授との共同研究である．本研究は，科学研究費補助金 基盤研究 (A)20H00576 研究代表者: 青嶋 誠「大規模複雑データの理論と方法論の革新的展開」，および，学術研究助成基金助成金 挑戦的研究 (萌芽) 22K19769 研究代表者: 青嶋 誠「テンソル構造をもつ巨大データの統計的圧縮技術の開発」から助成を受けたものである．
2010 Mathematics Subject Classification: 62H25, 62H30
キーワード：高次元統計解析, 非スパースモデリング∗e-mail: aoshima@math.tsukuba.ac.jp
web: https://www.math.tsukuba.ac.jp/˜aoshima-lab/jp/

が前提になっている．ところで，高次元データの母共分散行列に，スパース性は仮定できるのだろうか．

2.1. 高次元共分散行列の非スパース性
p 次非負定値対称行列 Σ を p 次元データの母共分散行列とし，その固有値を λ1 ≥ · · · ≥λp (≥ 0) とする．次の 3 つの高次元データセットについて，それぞれの固有値を上から10 個求めてみる．第 1 のデータセットは，Takeuchi et al. [11] による天体の分光マッピン
グによる波長データで，標本数は 231 地点，スペクトルは 1971 次元である．第 2 のデー
タセットは，von Roemeling et al. [12] による 2 群からなる淡明細胞型腎細胞癌 (ccRCC)のマイクロアレイデータで，ccRCC 群は 71 サンプル，健常群 (Normal) は 72 サンプルで，2 群とも遺伝子数 (プローブ数) が 54675 次元である．第 3 のデータセットは，Wanget al. [13] の次世代シーケンサーによる 2 群からなる加齢黄斑変性 (AMD) のゲノムデー
タで，AMD 群は 21 サンプル，健常群 (Normal) は 19 サンプルで，2 群ともゲノム領域
数が 3095656 次元である．図 1 は，寄与率 λj /tr(Σ) について，Yata and Aoshima [15] のノイズ掃き出し法による推定値 λ̃j /tr(S) をプロットしたものである 1 ．どの高次元デー
タも，最初の幾つかの固有値が飛び抜けて大きく，他の固有値は小さなものであることが見てとれる．特に，第 1 固有値は強くスパイクしており，一般に tr(Σ)/p → c (> 0)となることを考えれば，λ1 は O(p) のスピードで発散することが分かる．

図 1. 天体の分光マッピングによる波長データ (p = 1971)，2 つの群のマイクロアレイ
データ (p = 54675)，2 つの群の次世代シーケンサーによるゲノムデータ (p = 3095656)．
それぞれの高次元データについて，第 10 寄与率まで推定値をプロットした．
高次元データの母共分散行列の固有値に対して，Johnstone [9] や Paul [10] 等は，ランダム行列理論から標本固有値の漸近的性質を導出した 2 ．特に，スタンフォード大学の Johnstone 教授は，Tracy-Widom 分布に基づく標本固有値の極限分布を与えた ([9])．

1λ̃j は (4.2) 式で与えられるもの，S は標本共分散行列である．後述の (A-i) と (A-ii) のもと，λ̃j /tr(S) =λj {1 + oP (1)}/tr(Σ) (p → ∞, n → ∞) なる一致性をもつ．
2
これらの研究は，次元数 p と標本数 n が n/p → c (> 0) なる場合を考えている．高次元小標本データ(n/p → 0) には対応していないことに注意する．

そのとき，Σ の固有値には，次のようなスパイクモデルが仮定されていた．

λ1 , ..., λt は 1 よりも大きく，次元数 p には依存しない定数λt+1 = · · · = λp = 1(2.1)∑(2.1) のもとでは，∥Σ − I p ∥2F = ts=1 (λs − 1)2 = O(1) (p → ∞) となり 3 ，p × p 行列
Σ と p 次単位行列 I p の距離が有界になる．つまり，(2.1) は，Σ の非対角成分の殆どが零である，といった非常にスパースな共分散構造を想定していることになる．スパー
スモデリングでは，(2.1) と同様な「λ1 は次元数 p に関して有界である」という条件が仮定される．図 1 でも見たが，固有値が次元数 p に依存しないという仮定は，実データがもつ固有値構造から大きく乖離する．それでは，実データの固有値構造に飛び抜けて大きな固有値が現れる原因は何か．主たる原因は，高次元データを構成する変数間の相関にある．実データに対する Σ は，非対角成分の多くが非零となり，非スパースな共分散構造をもつのである．困ったことに，非スパースな共分散構造は高次元において巨大なノイズを生み，潜在情報を埋もれさせてしまう．これが，高次元データの解析を困難にしている．このような高次元データに，ランダム行列理論・スパースモ
デリング・深層学習などを，そのまま使用することは望ましくない．そこで，高次元データが非スパースな共分散構造をもつことを考慮して，新たに開発された方法論が非スパースモデリングである．

図 2. 高次元データの解析: スパース性と非スパース性．

2.2. 強スパイク固有値モデル
高次元データの母共分散行列 Σ は，次元数の増加とともに行列のサイズが大きくなり，
図 1 で見たように最初の幾つかの固有値は次元数に依存して大きなものとなる．特に，最大固有値 λ1 は，次元数の冪関数として捉えることが自然である．Aoshima and Yata
[4] は，高次元データの固有値モデルを 2 つに分類した．一つ目は，強スパイク固有値
モデル (strongly spiked eigenvalue (SSE) モデル) とよばれ，次のように定義される．
{ λ2 }
1
lim inf
>0
(2.2)
p→∞
tr(Σ2 )
いま一つは，弱スパイク固有値モデル (NSSE モデル) とよばれ，次のように定義される．

λ21→ 0 (p → ∞)tr(Σ2 )
3∥ · ∥F はフロベニウスノルムを表す．

(2.3)

任意の高次元データは，上記 2 つの固有値モデルの何れかに分類される．簡単のため，
以下のような固有値構造を考えてみる．

λs = cs pαs (s = 1, ..., t),λs = cs (s = t + 1, ..., p)(2.4)

ここで，cs > 0, α1 ≥ · · · ≥ αt > 0 は次元数 p に依存しない未知の実数，t は次元数 pに依存しない未知の自然数である．(2.4) は，α1 ≥ 0.5 のとき強スパイク固有値モデル(2.2) の一つとなり，α1 < 0.5 のとき弱スパイク固有値モデル (2.3) の一つとなる．
後述するが，弱スパイク固有値モデルにおいては種々の統計量に高次元漸近正規性という望ましい性質が得られるが，強スパイク固有値モデルでは高次元漸近正規性が成り立たない．それどころか，強くスパイクする固有値が一つでもあると，それが巨大なノイズとなって潜在情報を覆い隠してしまう．そのとき，強スパイク固有値モデ
ルは，非対角成分の多くが非零といった非スパースな Σ を想定したものになる．(2.2)
は母共分散行列に関する境界条件になっており，SSE モデルは非スパース性に対応し，NSSE モデルはスパース性に対応する．

3. 高次元漸近正規性
平均が p 次ベクトル µ，共分散行列が Σ の母集団から，n (≥ 3) 個のデータベクトル
x1 , ..., xn を無作為に抽出する．それらを並べて，p×n のデータ行列を X = [x1 , ..., xn ] とする．適当な直交行列 H = [h1 , ..., hp ] で Σ = HΛH T , Λ = diag(λ1 , ..., λp ) と分解する．
そのとき，X − [µ, ..., µ] = HΛ1/2 Z となる Z を Z = [z 1 , ..., z p ]T ，z s = (z1s , ..., zns )Tと表記する．Z の成分は 4 次モーメントが一様有界であると仮定し，さらに，必要な箇所で，適宜，次を仮定する．
2 2(A-i) E(zℓszℓt ) = 1，E(zℓs zℓt zℓs′ ) = 0，E(zℓs zℓt zℓs′ zℓt′ ) = 0，s ̸= t, s′ , t′(A-i) は，母集団の分布について，正規分布を緩めた仮定になっている．
母平均 µ に関する推測を考える．そのとき，共分散行列 Σ はノイズとなり，前述の∑SSE モデルと NSSE モデルはノイズモデルになる．標本平均ベクトル x = nℓ=1 xℓ /nについて，平均 2 乗誤差は E(∥x − µ∥2 ) = tr(Σ)/n となる．ここで，∥ · ∥ はユーク
リッドノルムである．一般に tr(Σ)/p → c (> 0) となることを考えれば，高次元小標本(p/n → ∞) においてE(∥x − µ∥2 ) → ∞となり，平均 2 乗誤差は発散してしまう．すなわち，次元数の増加とともに標本平均に含まれるノイズが巨大化し，x で µ を推定することは困難になる．ところで，∥x − µ∥2の分散は，(A-i) のもと p → ∞ のとき

(
)
p
2
)
tr(Σ2 )tr(Σ2 )tr(Σ2 )(n − 1) ∑ λ2s Var(zℓsVar(∥x − µ∥ ) = 2+=2 2 +O
n3
n3
n
n3
s=1
2

となる．さらに，球形条件 4 “tr(Σ2 )/{tr(Σ)}2 → 0 (p → ∞)”を満たすと仮定すると，
p → ∞ のとき (n は固定したままでも) 次が成り立つ．
√∥x − µ∥ = tr(Σ)/n{1 + oP (1)}4

tr(Σ2 )/{tr(Σ)}2 ≤ λ1 /tr(Σ) ≤ λ1 /{tr(Σ2 )}1/2 に注意すれば，NSSE モデル (2.3) は球形条件を満たす．

√
すなわち，高次元小標本 (p/n → ∞) において，x は中心 µ，半径 tr(Σ)/n の球面に収束する．図 3 は，単位行列 I p を共分散行列にもつ p 次元正規分布 Np (µ, I p ) について，
大きさ n = 3 の無作為標本による標本平均ベクトルを 200 回発生させ，x − µ を固有空
間 (ĥ1 , ĥ2 , ĥ3 で張られる空間) にプロットしたものである．ここで，ĥj は標本共分散行√√
列の第 j 固有ベクトルを表す．次元数が p = 4 のとき (左図) は，半径 p/n = 4/3 の球の周辺に点が散らばっている．一方，次元数が p = 1000 のとき (右図) は，ノイズが√√膨張して半径 p/n = 1000/3 の球の表面に点が集中し，球面集中現象が見てとれる．

図 3. 高次元小標本における x の球面集中現象．
球面での x の確率変動について，Aoshima and Yata [1, 3] は，次のような高次元漸近
正規性を証明した．
定理 1 ([3]). (A-i) を仮定する．NSSE モデルのもとで，p → ∞, n → ∞ のとき次が成り立つ 5 ．
∥x − µ∥2 − tr(Σ)/n L√−→ N (0, 1)
2tr(Σ2 )/n2L

ここで，−→ は分布収束を表す．

√
定理 1 は，高次元において，x が中心 µ，半径 tr(Σ)/n の球面の周りに正規分布で変動することを意味する．このことから，NSSE モデルであれば，球面付近の微小な確率変動を捉えることで統計的推測の精度を保証することが可能となる．しかしながら，
SSE モデルの場合，これが上手くいかない．強くスパイクする (非スパースな) ノイズ
が悪さをして，高次元漸近正規性は成立しない．いま，p = 1000 のとき，NSSE モデルとして Σ = diag(p1/3 , 1, ..., 1) を，SSE モデルとして Σ = diag(p, 1, ..., 1) を考えてみる．
図 4 は，p 次元正規分布 Np (µ, Σ) について，大きさの無作為標本による標本平√ n = 1022均ベクトルをもとに U = {∥x − µ∥ − tr(Σ)/n}/ 2tr(Σ )/n2 を計算し，これを 2000回発生させてヒストグラムを作成したものである．定理 1 の通り，NSSE モデルにおけるヒストグラムの形状は N (0, 1) に近いことが確認できる．一方，SSE モデルにおいては，高次元漸近正規性が成立しない．さらに，SSE モデルにおいては tr(Σ2 ) = O(λ21 )となるので，図 1 のように最大固有値が飛び抜けている場合，∥x − µ∥2 の分散が非常に大きくなる．つまり，SSE モデルが当てはまる状況では，高次元漸近正規性が壊れるだけでなく，推測の精度も著しく悪くなるのである．このことから，NSSE モデルを想定して (つまり，ノイズのスパース性を想定して) 構築された統計的推測法は，SSE
5

本講演で n → ∞ と表記するとき，n = log p 程度で十分である．つまり，高次元小標本であっても，結果は保証される．

図 4. NSSE モデル Σ = diag(p1/3 , 1, ..., 1) と SSE モデル Σ = diag(p, 1, ..., 1) の場合の，球面付近の x の漸近的挙動．U (p = 1000, n = 10) を 2000 回発生させ，ヒストグラムを作成した．実線は N (0, 1) の確率密度関数を表す．
モデルに対しては精度を保証することができず，それゆえ SSE モデルに使うわけにはいかない．それでは，SSE モデルも考慮した統計的推測法を構築することはできない
か．Aoshima and Yata [4] は，SSE モデルの場合は自動的に NSSE モデルに変換して推測を行うデータ変換法を考案した．次節でこれを説明するが，表記を簡単にするために，ここで SSE モデルを次のように簡略化しておく．

(A-ii) (次元数 p に依存しない) ある自然数 k に対して，(i) 1 ≤ s < s′ ≤ k のとき，lim inf (λs /λs′ − 1) > 0
p→∞

λ2
λ2
(ii) lim inf k > 0 かつ k+1 → 0, p → ∞p→∞ ΨkΨk+1(ただし，Ψj =∑p2
s=j λs である)

(2.4) の固有値構造であれば，k(< t) に対して αk ≥ 0.5 > αk+1 且つ 1 ≤ s < s′ ≤ k に対して cs ̸= cs′ ならば，条件 (A-ii) を満たす．(A-ii) は，図 1 の実データ解析をヒントに，
最初の k 個の固有値が飛び抜けて大きいことを表現した SSE モデルになっている．

4. データ変換法
Aoshima and Yata [4] は，次のようなデータ変換を考えた 6 ．
xℓ∗ = Axℓ (ℓ = 1, ..., n) (ただし，A =∑pT
s=k+1 hs hs (k ≥ 0) である)

(4.1)

∑p
T
変換後のデータについて，Var(xℓ∗ ) =s=k+1 λs hs hs (= Σ∗ とおく) となるので，
λmax (Σ∗ ) を Σ∗ の最大固有値とすれば，SSE モデル (A-ii) の条件 (ii) から
λ2k+1{λmax (Σ∗ )}2→ 0, p → ∞=Ψk+1
tr(Σ2∗ )が成立する．つまり，変換後のデータ xℓ∗ には，NSSE モデルが当てはまることになる．
∑k
データ変換に使う正射影行列は A = I p − s=1 hs hTs と書けるので，強くスパイクする
固有空間を如何に高精度に推定するかが鍵となる．本節では，高次元主成分分析 (PCA)として知られるノイズ掃き出し法を用いて固有値・固有ベクトルを推定する．
標本共分散行列S = (n − 1)−1 (X − X)(X − X)T
6(ただし，X = [x, ..., x] である)

強くスパイクする固有空間の個数 k は未知である．k の推定は，Aoshima and Yata [4] の S2.2 節を参照のこと．なお，k = 0 の場合は恒等変換となり，データは変換せずとも NSSE モデルが当てはまる．

について，S の固有値を λ̂1 ≥ · · · ≥ λ̂p (≥ 0)，各固有値 λ̂j に対する固有ベクトルを ĥj とする．ただし，ĥ1 , ..., ĥp は p 次元実ベクトル空間の正規直交基底をなすとする．S の固∑p
T
有値分解は S = s=1 λ̂s ĥs ĥs と書ける．S D = (n−1)−1 (X −X)T (X −X) とおくと，n次の対称行列 S D は S と正の固有値を共有する．S D の固有値を λ̂1 ≥ · · · ≥ λ̂n−1 (≥ 0)，
各 λ̂j に対する固有ベクトルを ûj = (û1j , ..., ûnj )T とする．ただし，û1 , ..., ûn−1 は互い∑n−1に正規直交とする．S D の固有値分解は S D = s=1 λ̂s ûs ûTs と書ける．S と S D の固有ベクトルには，次の関係がある．

ĥj = {(n − 1)λ̂j }−1/2 (X − X)ûj(j = 1, ..., n − 1)

Yata and Aoshima [14] は，固有値構造が (2.4) の場合に，λ̂j , ĥj , さらに主成分スコアについて，推定量の一致性 (と不一致性) を明らかにした．Aoshima and Yata [4] と Yataand Aoshima [16] は，一般の SSE モデルの場合にその理論を拡張している．SSE モデル(A-ii) を使って λ̂j と ĥj について簡単に述べると，次のようになる．
定理 2 ([4]). (A-i) を仮定する．SSE モデル (A-ii) のもとで，各 j (≤ k) について，p →∞, n → ∞ のとき次が成り立つ．
Tλ̂j /λj = 1 + δj + OP (n−1/2 ), (ĥj hj )2 = (1 + δj )−1 + OP (n−1/2 )∑p
ただし，δj = λ−1
js=k+1 λs /(n − 1) である．
上記の通り，λ̂j と ĥj は，δj というノイズをもつ．(2.4) の固有値構造であれば，αj = 0.5
Tに対して n = o(p0.5 ) だとすると δj → ∞ となり，λj /λ̂j = oP (1) や (ĥj hj )2 = oP (1) といった強不一致性を起こす．この問題を解決するために，Yata and Aoshima [15, 16] はノイズ掃き出し法という高次元における PCA を考案し，次のような固有値の推定量を与えた．
∑
tr(S D ) − js=1 λ̂s
λ̃j = λ̂j −(j = 1, ..., n − 2)
(4.2)
n−1−j
固有ベクトルの推定量は次のようになる．

h̃j = {(n − 1)λ̃j }−1/2 (X − X)ûj(j = 1, ..., n − 2)

(4.3)

Aoshima and Yata [4] と Yata and Aoshima [16] による次の定理は，δj → ∞ となる場合にも，ノイズ掃き出し法による推定量が一致性をもつことを示している．
定理 3 ([4]). (A-i) を仮定する．SSE モデル (A-ii) のもとで，各 j (≤ k) について，p →∞, n → ∞ のとき次が成り立つ．

λ̃j /λj = 1 + OP (n−1/2 ),T(h̃j hj )2 = 1 + OP (n−1 )データ変換 (4.1) を使うためには，Axℓ を構成する hTj xℓ (= xℓj とおく), j = 1, ..., kT

Tを推定する必要がある．Aoshima and Yata [4] は，ĥj xℓ や h̃j xℓ といった推定は非常に大きなバイアスが生むことを指摘し，次のような推定を考えた．
Tx̃ℓj = h̃ℓj xℓ(j = 1, ..., k)

(4.4)

ただし，
1/2

h̃ℓj = (n − 1)1/2 (X − X)ûℓj /{(n − 2)λ̃j },
ここで，

∑n

ûℓj = (û1j , ..., ûℓ−1j , −ûℓj /(n − 1), ûℓ+1j , ..., ûnj )Tℓ=1 h̃ℓj /n = h̃j となることに注意する．

5. 高次元 2 標本検定とデータ解析
母集団が 2 つあり，各母集団 πi (i = 1, 2) は平均に p 次ベクトル µi ，共分散行列に p 次∑p
非負定値対称行列 Σi (= s=1 λis his hTis ) をもつとする．次の検定を考える 7 ．

H0 : µ1 = µ2

vs. H1 : µ1 ̸= µ2(5.1)

各母集団から，ni (≥ 3) 個の p 次データベクトル xi1 , ..., xini を無作為に抽出し，データ
行列を X i = [xi1 , ..., xini ] とする．高次元小標本では標本共分散行列の逆行列が存在しないので，従来のホテリングの T 2 -統計量は使えない．高次元 2 標本検定では，次の統計量が使われる．

T = ∥x1 − x2 ∥ −22∑tr(S i )

ni

i=1

=2

2
∑
i=1

∑ni

T
ℓ<ℓ′ xiℓ xiℓ′ni (ni − 1)∑n1 ∑n2−2ℓ=1

Tℓ′ =1 x1ℓ x2ℓ′

n1 n2

(5.2)

∑ni
ただし，xi = ℓ=1
xiℓ /ni ，X i = [xi , ..., xi ], S i = (ni − 1)−1 (X i − X i )(X i − X i )T である．ここで，E(T ) = ∥µ1 − µ2 ∥2 となり, H0 のもとで
Var(T ) = 2

2∑

tr(Σ2i )4+
tr(Σ1 Σ2 ) (= K とおく)
n
(n
−
1)
n
n
i
i
1
2
i=1

となることに注意する．2 つの母集団がともに NSSE モデルのとき，T に関する高次元
漸近正規性が Chen and Qin [6] と Aoshima and Yata [1, 3] によって示され，それぞれ別の検定手法が与えられた．ところが，2 つの母集団のどちらか一方でも SSE モデルのとき，3 節の図 4 のように，T の高次元漸近正規性は成立しない．さらに，2 つの母集団がともに SSE モデルのとき，K = O(λ211 /n21 + λ221 /n22 ) となり，各母集団の最大固有値の影響で T の分散は非常に大きくなる．つまり，SSE モデルが当てはまる状況では，
T の高次元漸近正規性が壊れるだけでなく，検定の精度も著しく悪くなる．
Aoshima and Yata [4] は，SSE モデルも考慮し，データ変換法を使って高次元 2 標本検定を次のように与えた．各母集団のデータを，(4.1) と同様にデータ変換する．

xiℓ∗ = Ai xiℓ (ℓ = 1, ..., ni )(ただし，Ai =∑pT
s=ki +1 his his (ki ≥ 0) である)

変換後のデータを使うと，(5.2) 式の T は次のようになる．

T∗ =2
=2

ni
2 ∑
∑
xT Ai xiℓ′iℓ

∑n1 ∑n2−2ℓ=1

Tℓ′ =1 x1ℓ A1 A2 x2ℓ′

n (ni − 1)n1 n2
i=1 ℓ<ℓ′ i
∑ki
2 ∑ni
T
∑
ℓ<ℓ′ (xiℓ xiℓ′ −j=1 xiℓj xiℓ′ j )ni (ni − 1)∑n1 ∑n2
∑k1
∑k2
T
ℓ=1
ℓ′ =1 (x1ℓ −j=1 x1ℓj h1j ) (x2ℓ′ −j=1 x2ℓ′ j h2j )
−2
n1 n2
i=1

7

本節は，SSE モデルも考慮した高次元統計的推測を説明するために，基本的で本質的な高次元 2 標本検定を扱う．高次元判別分析などは，[5, 8] を参照のこと．なお，いわば基礎編となる NSSE モデルの高次元統計的推測については，青嶋・矢田 [2] が入門書となる．

ただし，xiℓj = hTij xiℓ である．ここで，E(T∗ ) = ∥A1 µ1 − A2 µ2 ∥2 (=∆∗ とおく) となり，
H0 のもとで2∑
4
tr(Σ2i∗ )+tr(Σ1∗ Σ2∗ )Var(T∗ ) = 2
ni (ni − 1) n1 n2
i=1(= K∗ とおく)∑pとなることに注意する．ただし，Σi∗ = s=ki +1 λis his hTis である．Aoshima and Yata [4]は，T∗ に関して高次元漸近正規性を証明した．さらに，Aoshima and Yata [4] は，T∗ を構成する hij と xiℓj を (4.3) 式の h̃ij と (4.4) 式の x̃iℓj で置き換えた T∗ の推定量 Tb∗ について，次の結果を与えた．
定理 4 ([4]). 2 つの母集団に (A-i) を仮定する．SSE モデルに (A-ii) を仮定する．さらに，
v = min{p, n1 , n2 } とおき，lim supv→∞ ∆2∗ /K∗ < ∞ を仮定する．適当な正則条件 8 のもと，v → ∞ のとき次が成り立つ．
√
Lb∗ −(Tb∗ − ∆∗ )/ K
→ N (0, 1)

b ∗ は Aoshima and Yata [4] の 5 節で与えられる K∗ の一致推定量である．
ただし，K
仮説 (5.1) は，次のように検定される．
√
bb ∗ > zα =⇒ 有意水準 α で H0 を棄却
T∗ / K(5.3)

ただし，zα は N (0, 1) の上側 α 点を表す．検定方式 (5.3) の第 1 種の過誤と検出力については，Aoshima and Yata [4] の定理 6 を参照のこと．
図 1 で紹介した Wang et al. [13] の次世代シーケンサーによる AMD ゲノムデータを解析する．ゲノム領域数が 3095656 (= p) 次元，AMD 群が 21(= n1 ) サンプル，健常群(Normal) が 19(= n2 ) サンプルであった．図 1 で見たように，両群とも SSE モデルが当てはまる．仮説 (5.1)√ を有意水準 α = 0.05 で検定する．検定方式 (5.3) の検定統計量をb ∗ = 23.44 > z0.05 = 1.65 となり，有意水準 5% で H0 は棄却され計算すると 9 ，Tb∗ / K
る．SSE モデル (非スパース性) を考慮し，データ変換法を使って非スパースなノイズを除去したことで潜在情報が浮き彫りとなり，2 群間の平均の差異が検出できた．
ちなみに，両群にスパース性を仮定してNSSE モデルのもとで検定統計量を計算する√
bと，T / K = 1.01 < z0.05 となり，H0 は棄却されない．しかし，この結論は統計的に保証されるものではない．なぜならば，NSSE モデル (スパース性の仮定) が妥当ではなく，検定統計量に高次元漸近正規性が成立していないからである．K = O(λ211 /n21 + λ221 /n22 )に注意すれば，T の潜在情報である E(T ) = ∥µ1 − µ2 ∥2 が非スパースなノイズとなる
最大固有値に押しつぶされ，2 群間の平均の差異を検出できなかったことが分かる．
次に，AMD 群と健常群 (Normal) の差異の原因となる変数 (ゲノム領域) を絞り込む．
いま，µ1 − µ2 = (µo1 , ..., µop )T ，x1 − x2 = (xo1 , ..., xop )T とおく．有意な変数の集合を
D = {j| |µoj | > 0} とおく．xo1 , ..., xop を |xo(1) | ≥ · · · ≥ |xo(p) | と絶対値の大きい順に並
b m = {j| |xoj | ≥ |xo(m) |} で D を推定する．適当な m とび替え，m 個の変数からなる D
8
9Aoshima and Yata [4] の定理 5 を参照のこと．
ki の推定については，[4] の S2.2 節を参照のこと．このデータセットは，(k̂1 , k̂2 ) = (7, 7) であった．

b m ) → 1 (p → ∞) となる 10 ．しかし，この変数選択
適当な正則条件のもと，P (D ⊆ Dはスパース性を当てにしたものである．
そこで，非スパース性を考慮し，データ変換法を使って xi∗ = Ai xi とする．(4.3) 式の h̃ij と (4.4) 式の x̃iℓj を用いて，xi∗ を
x̃i∗ = xi −ki∑

x̃ij h̃ij(ただし，x̃ij =∑niℓ=1 x̃iℓj /ni である)

j=1

¯o1 , ..., x̃¯op )T とし，x̃o1 , ..., x̃op を |x̃o(1) | ≥ · · · ≥ |x̃o(p) | と並びで推定する．x̃1∗ − x̃2∗ = (x̃e m = {j| |x̃oj | ≥ |x̃o(m) |} で D を推定する 11 ．
替え，m 個の変数からなる D
bm と De m を用いて，3095656 個の変数から m = 10, 100, 1000
結果を見てみよう．D個の変数をそれぞれ選択した．変数選択の精度を高次元 PCA で確認する．xij におい
b m の変数のみを利用した m 次元データを x
e m の変数のみを利用した m 次
b ij ，Dて, Dc =e ij とし，AMD 群と健常群 (Normal) を混合させたデータセットを X
元データを xf = (eb 1n1 , x
b 21 , ..., xb 2n2 )，X
e 1n1 , xe 21 , ..., xe 2n2 ) とおく．Yata and Aoshima(b
x11 , ..., xx11 , ..., x[17] は，規準化した主成分スコアを用いることで，高次元混合データのクラスタリングが可能であることを証明した．2 群間の差異に有意な変数を正しく選択できていれば，c
理論的には，主成分スコアを使って 2 群を精度良く分離できるはずである．そこで，X
f の各々について，n1 + n2 = 40 個のサンプルに対する（規準化した）第 1 主成分とX
em
スコアと第 2 主成分スコアを計算し，図 5 にプロットした．データ変換法を使った D
による変数選択は，2 群を精度よく分離できていることが見てとれる．データ変換法を使って非スパースなノイズを除去することで，3095656 個の変数から 2 群の差異に有意な 100 個の変数を絞り込めたのである．
非スパースモデリングの一幕をお見せした．データ解析に要した時間は，モバイル
PC で 5 分程度である．高い計算処理能力と大標本を前提とする深層学習が，経済的にも技術的にも環境的にも敷居の高いテクノロジーであるのに対して，高次元統計解析
による非スパースモデリングは，いわば庶民的な町工場の技術である．非スパースモ
デリングを使うことで，非スパースなノイズの影響を取り除くことができ，超高次元
データもたった数十程度の小標本で処理でき，モバイル PC でも高速かつ高精度に解析することができるのである．

参考文献[1] Aoshima, M. and Yata, K. (2011). Two-stage procedures for high-dimensional data, Sequential
Anal. (Editor’s special invited paper), 30, 356–399.
[2] 青嶋 誠，矢田和善 (2019). 高次元の統計学，共立出版，東京.
[3] Aoshima, M. and Yata, K. (2015). Asymptotic normality for inference on multisample, highdimensional mean vectors under mild conditions, Meth. Comput. Appl. Probab., 17, 419–439.
[4] Aoshima, M. and Yata, K. (2018). Two-sample tests for high-dimension, strongly spiked eigenvalue models, Statist. Sinica, 28, 43–62.
[5] Aoshima, M. and Yata, K. (2019). Distance-based classifier by data transformation for highdimension, strongly spiked eigenvalue models, Ann. Inst. Statist. Math., 71, 473–503.
[6] 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.
10

例えば，Aoshima and Yata [1] の 5 節を参照のこと．
b m よりも緩い条件で一致性をもつ．
Dm は D11 e

図 5. m = 10, 100, 1000 個の変数を選択した AMD 群と健常群 (Normal) の混合データに
b m : データ変換なし (左)，De m : データ変換法 (右)．
対する第 1 と第 2 の主成分スコア．D[7] Donoho, D.L. (2006). Compressed sensing, IEEE Trans. Information Theory, 52, 1289–1306.
[8] Ishii, A., Yata, K. and Aoshima, M. (2022). Geometric classifiers for high-dimensional noisy
data, J. Multivariate Anal. (Editor’s invited paper), 188, 104850.
[9] Johnstone, I.M. (2001). On the distribution of the largest eigenvalue in principal components
analysis, Ann. Statist., 29, 295–327.
[10] Paul, D. (2007). Asymptotics of sample eigenstructure for a large dimensional spiked covariancemodel, Statistica Sinica, 17, 1617–1642.
[11] Takeuchi, T, Yata, K., Aoshima, M., et al. (2022). High dimensional statistical analysis and its
application to ALMA map of NGC 253, arXiv preprint, arXiv:2203.04535.
[12] von Roemeling, C., Radisky, D., Marlow, L., et al. (2014). Neuronal pentraxin 2 supports clear
cell renal cell carcinoma by activating the AMPA-selective glutamate receptor-4, Cancer Res.,74, 4796–4810.
[13] Wang, J., Zibetti, C., Shang, P., et al. (2018). ATAC-seq analysis reveals a widespread decreaseof chromatin accessibility in age-related macular degeneration, Nat. Commun., 9, 1364.
[14] 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.
[15] Yata, K. and Aoshima, M. (2012). Effective PCA for high-dimension, low-sample-size data withnoise reduction via geometric representations, J. Multivariate Anal., 105, 193–215.
[16] Yata, K. and Aoshima, M. (2013). PCA consistency for the power spiked model in highdimensional settings, J. Multivariate Anal., 122, 334-354.
[17] Yata, K. and Aoshima, M. (2020). Geometric consistency of principal component scores forhigh-dimensional mixture models and its application, Scand. J. Stat., 47, 899–921.