2015 年度科学研究費シンポジウム

大規模複雑データの理論と方法論：最前線の動向
科学研究費補助金 基盤研究 (A) 15H01678「大規模複雑データの理論と方法論の総合的研究（研究代表者：青嶋誠）」，学術研究助成基金助成金 挑戦的萌芽研究 26540010「ビッグデータの統計学: 理論の開拓と 3V への挑戦（研究代表者：青嶋誠）」によるシンポジウムを下記のように催しますので，ご案内申し上げます．
誠(筑波大学)

矢田和善(筑波大学)

日野英逸(筑波大学)

青嶋記
日時：2015 年 11 月 16 日 (月)〜18 日 (水)場所：筑波大学自然系学系棟Ｄ棟 D509 (筑波キャンパス内)〒305-8571 茨城県つくば市天王台 1-1-1

プログラム
11 月 16 日 (月)14:00〜14:05
14:05〜14:35

開会
園田

翔

（早稲田大・先進理工学研究科）

村田昇（早稲田大・先進理工学研究科）

ReLU ネットワークの積分表現理論
14:45〜15:25

矢田

和善 （筑波大・数理物質系）

青嶋誠（筑波大・数理物質系）

High-dimensional two-sample tests in general settings
15:40〜16:10

江口

翔一 （九州大・数理学府）

従属型一般化線形モデルにおけるモデル評価

16:20〜17:00

福地純一郎（学習院大・経済学部）

Nonparametric bootstrap for Gupta’s selection procedure
11 月 17 日 (火)9:20〜9:50

中川智之 （広島大・理学研究科）

高次元の場合での誤判別確率の推定について

10:00〜10:30

岡

紘之

白石

博

（慶應義塾大・理工学研究科）（慶應義塾大・理工学部）

高次元の下での有効フロンティアの統計的推定

10:45〜11:25

蛭川

潤一 （新潟大・自然科学系）

坂井俊介 （新潟大・自然科学研究科）

Rank tests for an ARMA model against other tv-ARMA models
11:35〜12:15

星野伸明 （金沢大・経済学類）

疎な分割表上の分布

12:15〜13:45

昼食

13:45〜14:25

柿沢佳秀 （北海道大・経済学研究科）

Some multivariate asymmetric kernel estimators for density functions
14:35〜15:15

柳本

武美 （統計数理研究所）

事前密度のベイズ尤度とその含意

15:30〜16:10

若木宏文 （広島大・理学研究科）

ランダム係数を持つ GMANOVA モデルの変数選択規準

16:20〜17:00

筑瀬靖子 （香川大）Peter E. Jupp（University of St. Andrews）
Distribution theory and its asymptotics in shape analysis
18:00〜

懇親会

11 月 18 日 (水)9:10〜9:40

茶之原

直人（広島大・理学研究科）

中川智之 （広島大・理学研究科）

若木宏文 （広島大・理学研究科）

修正コレスキー分解を用いた分散共分散行列の縮小推定について

9:50〜10:20

渋谷

明（東京理科大・理工学研究科）

生亀

清貴 （東京理科大・理工学部）

富澤

貞男 （東京理科大・理工学部）

順序カテゴリ正方分割表における対角指数対称モデルについて

10:35〜11:05

稲津佑（広島大・理学研究科）

パラレルプロファイルモデルにおけるランダムエフェクト共分散構造に対する尤度比統計量の高次元漸近展開

11:15〜11:55

鎌谷

研吾 （大阪大・基礎工学研究科, JST CREST）

マルコフ連鎖モンテカルロ法の高次元漸近論

11:55〜12:00

閉会報告書：ReLU ネットワークの積分表現理論 ∗†
園田
早稲田大学

翔‡

村田

昇

先進理工学研究科

Sonoda Sho Noboru Murata
Faculty of Science and Engineering
Waseda University
本発表では，まずディープラーニングの現状について簡単に説明したのち，ReLU ネッ
トワークを具体例として積分表現理論について紹介した。
ディープラーニングでは，従来よりも多くの隠れ層を備えたニューラルネットを学習させる。隠れ層を重ねることで，情報表現は階層化され，情報の流れは組合せ的に複雑になる。この原理によって，浅いニューラルネットよりも効率的に学習ができることは，以前から予見されていた。しかし現実には，隠れ層が深くなるにつれて誤差信号が指数的に減衰するので，バックプロパゲーションは実行不可能であった。ディープニューラルネット
では，モデルパラメータの数は数百万から数億のオーダーにのぼることもあるのに対し，
データサイズは多くとも数千から数万のオーダーに留まる。それにもかかわらず，画像認
識を中心としとする応用分野で好成績を収めていることは，一種の驚きである。
ニューラルネットの積分表現とは，ニューラルネットを連続化した関数∫
g(x) =
Rm ×R

T(a, b) η(a · x − b)dadb,

x ∈ Rm(1)である。適当な方法で積分表現を離散化することで，通常のニューラルネットを得る。積分表現 (1) は，T(a, b) の η による双対リッジレット変換として知られる。T として適当な関数 f のリッジレット変換∫Rψ f (a, b) :=Rmf (x) ψ(a · x − b)dx,∗ 「大規模複雑データの理論と方法論：最前線の動向」† 本研究は JSPS 科研費 15J07517 の助成を受けたものです。
‡ s.sonoda0110@toki.waseda.jp

1(a, b) ∈ Rm × R(2)

をとると，f, ψ, η に対する適当な条件のもとで再生公式

Rη† Rψ f (x) = f (x),(3)

が成り立つ。再生公式は，ニューラルネットが万能であることを表している。
積分表現を用いて，活性化関数が ReLU のように非有界な場合にもニューラルネットの万能性が保たれることを説明した。形式的には，積分表現が b に関して畳込み積分になっていることを利用して，畳込みに対する微分の性質 (部分積分)

∫

Rm ×R

∫

′

T(a, b)η (a · x − b)dadb =Rm ×R∂b T(a, b)η(a · x − b)dadb,

(4)

から従う。
再生公式から，ニューラルネットが学習の結果として獲得する情報表現の正体がリッジ
レット変換であることを説明した。隠れ層における情報表現 T = Rψ f は，本質的に一意性がない。なぜならば，ひとつの与えられた η に対して許容条件を満たす ψ を複数構成できるためである。つまり隠れ層における情報表現 T = Rψ f には，ψ の不定性がある。
可能な（許容条件を満たす）ψ からさらに「よい」ψ を選択する方法として，積分表現

Rη† Rψ f において Rψ f が係数に相当することに基いて，∥Rψ f ∥1 を最小化することでス
パーシティの高い係数を選ぶ方法を紹介した。

min ∥Rψ f ∥1 .

ψ∈Aη

(5)

通常のニューラルネットとして考えると，この定式化は L1 正則化付きの最小二乗法に相当する。ただし積分表現の式では，常に許容条件を満たしていることが前提なので，二乗誤差は適当な近似誤差 ε > 0 で押さえられたうえでの最適化として考えてよい。

min ∥f −aj ,bj ,cj∑

cj · η(aj · x − bj )∥2 + λ∥c∥.

(6)j

積分表現 Rη† Rψ (a, b) を直接数値計算することで，バックプロパゲーションを経ずに学習する方法を紹介した。本発表では，ランダムに離散化を行うオラクル分布という方法と，数値積分を計算する方法を説明した。
ディープニューラルネットに対する積分表現はまだない。本発表では，例えば隠れ層が2 層の場合でも，以下のように入れ子構造になるために難しいということを説明した。
∫
g(x) =T′ (a′ , b)η(a′ · h(x) − b)da′ db,
(7)
′
Rm ×R
∫
h(x) =
T(a, b)η(a · x − b)dadb.
(8)
Rm ×R

2

High-dimensional two-sample tests in general settings
矢田 和善
青嶋 誠

(筑波大数理物質)(筑波大数理物質)

1. はじめにゲノム科学等に見られる高次元データの一つの特徴は，データがもつ次元数の膨大さにある．次元数は標本数を遥かに超える．いま一つの特徴は，高次元データは豊富な情報を有するものの，それが巨大なノイズに埋もれて見つけ難いことである．これ
らの理由から，通常の多変量解析法では高次元データの推測に精度を保証することができず，間違った解析結果を導くことさえある．
高次元二標本問題について，Bai and Saranadasa (1996, Statist. Sinica)，Chen andQin (2010, Ann. Statist.) や Aoshima and Yata [1] 等は，母集団分布が正規性よりも緩い仮定のもとで，ユークリッド距離に基づく統計量を与え，その漸近正規性を示し，それに基づく二標本検定法を考えた．一方で，マハラノビス距離に基づく統計量については，Cai et al. (2014, JRSS-B) が高次元の枠組みで扱っている．
本講演では，二つの固有値モデル，strongly spiked eigenvalue (SSE) モデルと nonSSE (NSSE) モデルに基づいて高次元二標本問題を考えた．まず，固有値の発散速度が遅い NSSE モデルのもと，高次元検定統計量の漸近一致性と漸近正規性を与えた．それらに基づき，ある種の最適性について議論し，ホテリングの T 2 検定等のマハラノビス
距離に基づく統計量が最適性を持たないことを示した．さらに，検定方式を選ぶための理論的な指標も与えた．次に，固有値の発散速度が速い SSE モデルのもと，新たな高次元検定統計量を提案し，その漸近正規性を示した．それらに基づく検定法が，SSE
モデルのもと他の検定法より優れていることを理論的かつ数値的に示した．

2. 高次元検定統計量に関する漸近的性質
母集団が 2 個あると想定し，各母集団 (πi ) は平均に p 次のベクトル µi ，共分散行列に p 次の正定値対称行列 Σi (> O) をもつと仮定する．いま，各母集団 πi から ni (≥ 4)
個の p 次元データ xi1 , ..., xini を無作為に抽出する．母集団 πi , i = 1, 2 の高次元データに，xij = Γi wij + µi (j = 1, ..., ni ) なるモデルを考える．ここで，Γi は Γi ΓTi = Σi な
る p × ri 行列，wij は E(w ij ) = 0, Var(wij ) = I ri とし，wij = (wi1j , ..., wiri j )T の各成
分は 4 次モーメントが一様有界と仮定する．ただし，I ri は ri 次元の単位行列である．
母集団 πi の分布には，必要な箇所でその都度，次を仮定する．

(A-i)

2
2
) = 1，E(wiqj wisj witj wiuj ) = 0，q ̸= s, t, u; i = 1, 2.
wisj
E(wiqj(A-i) は正規分布を緩めた仮定になっている．高次元のもと，
H0 : µ1 = µ2

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

なる高次元平均ベクトルの検定を考え，漸近的な最適性を議論する．
いま，非負定値対称行列 A を考える．各 i で，標本平均を xini ，標本共分散行列を
S ini とし，(1) の検定統計量をT (A) = (x1n1 − x2n2 ) A(x1n1 − x2n2 ) −
T

2
∑
i=1

tr(S ini A)/ni

とおく．T (I p ) は，Aoshima and Yata [1] 等で与えた，ユークリッド距離に基づく統計量と一致することに注意する．ここで，µA = A1/2 (µ1 − µ2 )，Σi,A = A1/2 Σi A1/2
とおく．そのとき，E(T (A)) = ||µA ||2 (= ∆(A) とおく)，Var(T (A)) = K(A) +∑4 2i=1 µTA Σi,A µA /ni となる．ただし，K(A) = 22∑tr(Σ2i,A )
i=1

ni (ni − 1)+4tr(Σ1,A Σ2,A )n1 n2である．いま，任意の正定値対称行列 M の最大固有値を，λmax (M ) で表すとし，Σi,Aには次の NSSE モデルを仮定する．

λmax (Σi,A )→ 0, p → ∞; i = 1, 2.
tr(Σ2i,A )1/2

(A-ii)

さらに，m = min{p, n1 , n2 } とおき，次のどちらかを仮定する．

(A-iii)K(A)→ 0, m → ∞;{∆(A)}2(A-iv) lim supm→∞{∆(A)}2< ∞.
K(A)H0 のもと，(A-iv) を満たすことに注意する．このとき，次の定理を得た．
定理 1.

(A-iii) を仮定する．m → ∞ のとき次が成り立つ．
T (A)= 1 + oP (1).
∆(A)

定理 2. (A-i) を仮定する．さらに，(A-ii) と (A-iv) のもとで，m → ∞ のとき次が成り立つ．
T (A) − ∆(A)√
⇒ N (0, 1),
K(A)ここで，⇒ は分布収束を表す．
定理 1 と 2 に基づき，各種 A による検定法の高次元における特徴と性質を報告し，NSSE モデルのもと漸近的な最適性を議論した．それにより，A が共分散行列の逆行列
であるようなマハラノビス距離に基づく統計量が最適性を持たないことを示し，検定方式を選ぶための理論的な指標を与えた．さらに，固有値の発散速度が速い SSE モデルを与え，新たな高次元検定統計量を提案した．SSE モデルのもと，その漸近正規性を与えることで検定法を構築し，それが他の検定法より優れていることを理論的かつ数値的に示した．

参考文献[1] Aoshima, M. and Yata, K. (2015a). Asymptotic normality for inference on multisample,
high-dimensional mean vectors under mild conditions. Methodol. Comput. Appl. Probab.
17, 419-439.
[2] Aoshima, M. and Yata, K. (2015b). Two-sample tests for high-dimension, strongly spiked
eigenvalue models, submitted.

従属型一般化線形モデルにおけるモデル評価
九州大学大学院数理学府

江口 翔一

Bayes 因子に基づいた相対モデル記述評価においては, データの対数周辺尤度の漸近挙動に基づいて最適モデルを決定する. すなわち, 対数周辺尤度の主要項を使いやす
い形で近似することができれば記述的モデル評価規準として用いることができる. このとき, 適当な正則条件および Taylor 展開により古典的な Bayes 情報量規準 BIC が導出されるが, 従属データモデルなどでは, BIC 導出の数学的正当性が厳密に議論がされない場合がしばしばある. そこで, 一般の疑似対数尤度関数を対象として, 漸近混合
正規性やモデル誤特定の場合まで視野に入れたモデル設定で BIC の拡張を行うことは理論上重要である. 本シンポジウムでは特に, 一種の従属型一般化線形モデルに対して, Schwarz 型統計量によるモデル評価の正当性を与えた. 本結果により, Lv and Liu(2014, Section 3.1) の結果が大幅に拡張された.
′′

Yn = (Y1 , . . . , Yn ) を n 次元確率変数, Xj = (xj,1 , . . . , xj,p ) (j = 1, . . . , n),′

Xn = (X1 , . . . , Xn ) を n × p random time series とする. ここで, (Y1 , ..., Yn ) が(X1 , ..., Xn )-条件付き独立であり, 各 Yj の (X1 , ..., Xn )-条件付き分布は Xj にのみ依存し, Xn の分布については未知でよいが θ には依らないものとする. このとき, モ
デル選択の対象は (Xn , Yn ) の同時分布であるが, 候補モデルにおける確率 (密度) 関
数は

f (xn , yn ; θ) = f (xn )f (yn |xn ; θ) = f (xn )n∏f (yj |xj ; θ)j=1と表すことができ, 条件付き分布 Yn |Xn についてのみ考察すればよい. この条件付き分布に一般化線形モデルを仮定すると, 疑似対数尤度関数はn∑()
Hn (θ) =Yj Xj′ θ − b(Xj′ θ)
j=1を考えればよい. ここで, b(·) は凸関数である (e.g. 線形 Gauss 回帰 : b(θ) = |θ|2 /2,
線形 Logistic 回帰 : b(θ) = log(1 + eθ )). この疑似対数尤度関数に対して対数周辺疑似尤度の確率展開を導出するための条件を与え, 特に, Xn については広範な時系列モデルを扱える一般的な条件を与えた. このときパラメータ θ ∈ Θ (⊂ Rp ) の事前分布

を π とし, 適当な仮定の下, 対数周辺疑似尤度の確率展開が導出できた:

(∫
log

)
exp{Hn (θ)}π(θ)dθΘpp
log n + log 2π + log π(θ̂n )
2
2
)
( ∑
n
1
1
2
′
′
∂ b(Xj θ̂n )Xj Xj + op (1)− log det2
n j=1

= Hn (θ̂n ) −ただし, θ̂n は疑似最尤推定量, つまり θ̂n ∈ argmaxθ∈Θ Hn (θ) である. この発散項に着目することにより, Bayes モデル (π, Hn ) に対する疑似 Bayes 型情報量規準(Quasi-BIC; QBIC) が
QBIC := −2
= −2

n
∑
(

j=1
n
∑

Yj Xj′ θ̂n − b(Xj′ θ̂n ))+ log det

(∑
n

∂

2

b(Xj′ θ̂n )Xj Xj′

)

j=1

(
)
Yj Xj′ θ̂n − b(Xj′ θ̂n ) + p log n + Op (1)j=1という形で得られた. これにより, Schwarz 型統計量 QBIC (BIC) の使用を正当化できた. 候補モデルごとに QBIC を計算し, それが最小となるモデルを最適なモデルとして選択する.
この QBIC と BIC, 赤池情報量規準 AIC を用いてモデル選択の実験を行い, その結果を比較した. すると, QBIC は BIC と AIC の中間となるような振る舞いをみせ,
データの増加とともに最適なモデルを選択する回数は増加する傾向がみられた.
参考文献[1] Lv, J. and Liu, J. S. (2014), Model selection principles in misspecified models.
J. R. Stat. Soc. Ser B. Stat. Methodol., 76, 141–167.

Nonparametric Bootstrap for Gupta’s Selection Procedure
福地純一郎

学習院大学

2015 年 11 月

1

母集団選択問題

Π1 , . . . , Πk は k 個の母集団で, 平均がそれぞれ µ1 , . . . , µk , 共通の分散 σ 2 を持ち, 分布関数はそ
れぞれ F (x − µ1 ), . . . , F (x − µk ) であるとする. ただし F は未知であるとする. 母集団 Πi (i =1, . . . , k) からの標本を {Xij : j = 1, . . . , n} で表すとする. X i を標本 {Xij , j = 1, 2, . . . , n} の標本
平均, s2n をプールされた標本分散とする. 母平均のパラメータ空間を Ω = {(µ1 , . . . , µk ) : µi ∈ R}
として, 順序付けた母平均を µ[1] ≤ µ[2] ≤ · · · ≤ µ[k] と表す. また, X (i) は母平均が µ[i] である母集団からの標本の標本平均を表すとする. また順序付けた標本平均を X [1] ≤ X [2] ≤ · · · ≤ X [k]と表す. Gupta (1965) が正規母集団の場合に提案した母集団選択ルールは以下で定義される.
母集団選択ルール R :dsn
X i ≥ X [k] − √n(1)となる母集団 Πi を選択する.
記号 CS で正しい選択, つまり最良の母集団が選択される部分集合に含まれる事象を意味す
るものとする. このとき正しい選択の確率 P (CS|R) の Ω 上の下限が µ1 = · · · = µk のときに与
えられ,
(√
)
√
inf P (CS|R) =Pn(X (j) − µ[j] ) ≤ n(X (k) − µ[k] ) + dsn , j = 1, . . . , k − 1(2)Ωであることが示される. Gupta の部分集合選択の方法では, 定数 P ∗ , (1/k < P ∗ < 1) に対して
inf Ω P (CS|R) = P ∗ となる d を用いる. 定数 d = d(P ∗ , n, k, F ) は確率変数
(
)
√
√
′
′
′−1
Wn = snmax
n Xi − n Xk .
1≤i≤k−1∑n′iidの P ∗ 分位点であることがわかる. ただし Xij′ ∼ F (x), X i = n−1 j=1 Xij′ , (i = 1, . . . , k),)∑k ∑n ( ′
′ 2
−1
−
X
=
((n
−
1)k)
s′2
X
. である. 分布関数 F は未知であるから d を推定す
ij
n
i
j=1
i=1
る必要がある.

2

ブートストラップ法

ブートストラップ法を用いて d を推定する以下のような方法の性質を調べた. 1. Fn を中心化標
本 Xn = {Xij − X i : i = 1, . . . , k, j = 1, . . . , n} の経験分布関数とする. 2. Xn を与えたときの条件付き Fn からの無作為標本 {Xij∗ : i = 1, . . . , k, j = 1, . . . , n} をブートストラップ標本とする. 3. ブートストラップ平均, ブートストラップ分散, Wn のブートストラップ・バージョン Wn∗を求める. 4. 2, 3 の手続きを多数 (B) 回繰り返すことにより, Wn∗ の分布の P ∗ 分位点を求める.
この方法について, 以下の結果が成り立つ.

1

iid

Theorem 1 Xij ∼ F (x − µi ), (i = 1, . . . , k, j = 1, . . . , n) であり, E(Xij ) = µi , V(Xij ) = σ 2を仮定する. このとき supx∈R |P (Wn∗ ≤ x|Xn ) − P (Wn ≤ x)| → 0, w.p.1. が成り立つ.
次に Wn の分母を σ で置き換えた W̃n の分布に関して以下が成り立つ.

Theorem 2 確率変数列 X1 , . . . , Xn は i.i.d. 確率変数列で非格子点分布を持つとする
. また,
√
2
3
−1
E(X1 ) = µ, V(X1 ) = σ , E(|X1 | ) < ∞ を仮定する. Gn は σn(X − µ) の確率分布関数,
33
γ = E(X1 − µ) /σ , k は正の整数,
Qk (x) = 6−1 γkΦk−1 (x)(1 − x2 )ϕ(x), x ∈ R.とする. このとき x について一様に Gkn (x) = Φk (x) + n−1/2 Qk (x) + o(n−1/2 ) が成り立つ.

Theorem 3 確率変数列 Xi1 , . . . , Xin , (i = 1, . . . , k) は i.i.d. 確率変数列で非格子点分布を持つとする. また, E(X11 ) = µ, V(X11 ) = σ 2 , E(|X11 |3 ) < ∞ を仮定する. また E(X11 ) = 0,
V(X11 ) = σ 2 , E(|X11 |3 ) < ∞. W̃n を
)
(
√
√
1
n Xi − n Xk .
W̃n =
max
σ 1≤i≤k−1() ∫∞と定義する. このとき P W̃n ≤ x = −∞ Φk−1 (y + x) dΦ(y) + n−1/2 Jn (x; γ) + o(n−1/2 ) が x
について一様に成り立つ. ただし
∫ ∞
(
)
}
{
−1
Jn (x; γ) = 6 γ(k − 1) 1 − (y + x)2 ϕ(y + x) + (y 3 − 3y)Φ (y + x) Φk−2 (y + x)ϕ(y)dy.
−∞

である.

Theorem 4 確率変数列 Xi1 , . . . , Xin , (i = 1, . . . , k) は i.i.d. 確率変数列で非格子点分布を持つとする. また, E(X11 ) = µ, V(X11 ) = σ 2 , E(|X11 |3 ) < ∞ を仮定する. また E(X11 ) = 0,
V(X11 ) = σ 2 , E(|X11 |3 ) < ∞. F の特性関数 F̂ は Cramér 条件
lim F̂ (t) < 1.

∥t∥→∞満たすと仮定する. このとき
(
) ∫ ∞
∗
P W̃n ≤ x Xn =
Φk−1 (y + x) dΦ(y)+ n−1/2 Jn (x; γ̂n ) + o(n−1/2 ),

(3)

−∞

∑k ∑n
が確率 1 で x について一様に成り立つ. ただし γ̂n = (nk)−1 i=1 j=1 (Xij − X i )3 /σ̂n3 . である.
(3) 式から()()
sup P W̃n∗ ≤ x Xn − P W̃n ≤ x = o(n−1/2 ), w.p.1.
x∈Rが成り立つことがわかる.

References[1] Babu, G.J. and Singh, K (1984). On one-term Edgeworth correction by Efron’s boostrap
Sankya Ser. A, 46, 219-232[2] Bickel, P.J. and Freedman, D.A. (1981). Some asymptotic theory for the bootstrap. Ann.
Statist. 9, 1196-1217.
[3] Gupta, S.S. (1965). On Some Multiple Decision (Selection and Ranking) Rules. Technometrics, 7, No. 2, 225-245[4] Singh, K. (1981). On asymptotic accuracy of Efron’s bootstrap. Ann. Statist. 9, 1187-1195.
2高次元の場合での誤判別確率の推定について
広島大学大学院

理学研究科

中川 智之

本報告は誤判別確率の推定について行った. 誤判別確率の推定方法としては Cross Validation(CV, Lachenbruch and Mickey(1968), Stone(1974)) や漸近展開 (Okamoto(1963, 1968),
McLachlan(1974), Fujikoshi and Seo(1998), etc.) などがある.
本研究は 2 母集団 Π1 , Π2 の判別分析を考える. Πi の個体を誤って Πj に判別してしまう確率 (誤判別確率) を P (i|j) と表す. p 次元母集団 Πk からのトレーニングデータ
Xk1 , ..., XkNk (k = 1, 2) とする. このとき, 線形判別関数を
{
(
) }
(
)′
L(X) = L(X; X̄1 , X̄2 , S) = X̄1 − X̄2 S −1 X − X̄1 + X̄2 /2 ,とする. ここで, X̄k , S は次で定義される不偏推定量である. 線形判別関数に関する漸近展開としては, Okamoto(1963, 1968) で大標本漸近理論 (Nk → ∞, p:fix) の場合での正規母集団に対する漸近展開が次の形で与えられている.

P (2|1) = P (L(X) ≤ 0|X ∈ Π1 ) = Φ (−∆/2) +

1
1
1
P1 +
P2 +
P3 + O2
N1
N2
N −2

(1)

ここで, Om は N1−1 , N2−1 , p−1 に関するオーダー, ∆2 = (µ1 −µ2 )′ Σ−1 (µ1 −µ2 ), N = N1 +N2である. また, McLachlan(1974) では漸近展開式に推定量を入れたときのバイアスの補正を行っている. Deev(1970) では N1 = N2 の場合の漸近近似式, Fujikoshi and Seo(1998) では
N1 ̸= N2 の場合を含む漸近近似式がそれぞれ高次元大標本漸近理論 (Nk , p → ∞, p/N →c0 ∈ (0, 1)) の場合で与えられている. また, Tonda and Wakaki(2003, TR) では漸近展開
が次の形で与えられている.

P (2|1) = Φ (ν) + ϕ(ν)f1 (∆) + O3/2 ,
)
(
p p
,
ν = ν ∆,
N1 N2(
)1/2 {}{}−1/2
1 N −p
p
pN
2
2
=−
∆ +
(N1 − N2 )∆ +2 N −1
N1 N2
N1 N2

(2)

また, 推定量を入れたときのバイアス補正もしている. 一方で, CV での推定に関しては大標本漸近理論の場合で, バイアスが O2 であることが知られている.
本報告では高次元大標本漸近理論の場合での CV に関する漸近的評価とバイアス補正の提案を行った. この章では線形判別関数での判別の誤判別確率を leave-one-out CV で推定したときのバイアスと MSE の評価を行う. まず,(j)X̄1 , S (j) : X1j を除いた推定量,(j)

Lj (X1j ) = L(X1j , X̄1 , X̄2 , S (j) ),とおくと, leave-one-out CV の推定量は,
11 ∑(2|1) =
1(Lj (X1j ) ≤ 0)
N1 j=1NP̂

CVと書ける. ここで, 1(·) は定義関数である. バイアスは (2) より, 次が成り立つ.
[]
E P̂ CV (2|1) − P (2|1) = O1 ,大標本漸近理論の下ではバイアスは O2 であったが, 高次元大標本漸近理論の場合ではオー
ダーが 1 つ上がることがわかる. 次に MSE の評価を行う. MSE は分散とバイアスの 2 乗に分解できる. すなわち,
[{
}2 ]
(
)
[
]2
CV
CV
CV
E P̂ (2|1) − P (2|1)= Var P̂ (2|1) + E P̂ (2|1) ,なので, P̂ CV (2|1) は漸近不偏推定量であり, 正規分布の仮定から()
CV
Var P̂ (2|1) → 0
が成り立つことが示せる. したがって, 以下が成り立つ.
[{
}2 ]
CV
E P̂ (2|1) − P (2|1)→ 0, (Ni → ∞, p → ∞, p/N → c0 ∈ (0, 1)) .
また, 本報告では k-hold CV で推定した推定量 P̂ CVk のバイアスは O(k −1 ) であるので,Yanagihara and Fujisawa(2012) での考え方を応用して, k-hold CV に関してバイアス補正の方法を提案した. まず, N1 = km となると仮定する. P̂new (2|1) 以下のように定義する.
{)}
k1 − k2 ( CVk
CVk1
CVP̂ 2 (2|1) − P̂ k1 (2|1)
P̂new (2|1) = P̂(2|1) −
k1 − 1
そのとき, 次のオーダーまでバイアス補正ができる.
[]
E P̂new (2|1) − P (2|1) = O(k1−2 )また, 最後に数値比較を行った.

参考文献[1] Fujikoshi, Y. , Seo, T. (1998), Asymptotic approximations of EPMC’s of the linear
and the quadratic discriminant functions when the sample sizes and the dimension
are large, Random Oper. Stochastic Equations, 6, 269-280.
[2] Lachenbruch, P. A. and Mickey, M. R. (1968), Estimation of Error Rates in Discriminant analysis. Technometrics, 10, 1-11.
[3] Stone, M. (1974), Cross-Validatory choice and assessment of statistical predictions.
J. R. Statist. Soc. , B36, 111-147.
[4] Tonda, T. and Wakaki, H. (2003), EPMC Estimation in Discriminant Analysis When
the Dimension and Sample Sizes are Large. Hiroshima Statistical Research Group
Technical Report, TR, 03-08.
[5] Yanagihara, H. and Fujisawa, H. (2012), Iterative bias correction of the crossvalidation criterion. Scand. J. Stat. , Vol. 39, 116-130.

高次元の下での有効フロンティアの統計的推定
慶應義塾大学理工学研究科 岡 紘之慶應義塾大学理工学部 白石 博

1 研究背景ポートフォリオ理論は，リスク回避的な投資家が分散投資を行い，自身のポートフォリオにおいてポートフォリオ収益率を高めるためには，どのように最適化すればよいかを決定するための理論である．実際には，資産のリターンはランダムな変数であり，どのような平均と分散をもつ分布に従っているかは未知であるため，サンプルから推定することを考える．近年の株式市場などでは，市場の規模が大きくなっており，分散投資する資産の数が膨大となっている．一般に，d が大きくなるにつれて大標本漸近理論での近似が悪くなってしまうことが知られており，このような状況を鑑みて，高次元データにおける解析を考える必要性が主張されている．
本研究では，ポートフォリオによって実現化される点の領域（実現可能領域）の左側境界を表す有効フロンティア（eﬃcient frontier）を次元 d，標本サイズ n のデータ行列 X = (X1 , . . . , Xn ) から推定し，

d → ∞,

n → ∞,d→ ρ ∈ (0, 1)
n

(1.1)

という設定のもとでの漸近的挙動について結果を報告した．

2 推定問題X1 , . . . , Xn ∈ Rd を平均ベクトル µ，分散共分散行列 Σ をもつある分布からの無作為標本とする．ただし，Σ は正定値行列であるとし，Xj は 4 次までのモーメントを持つと仮定する．また，標本平均 X̄n と標本分散共分散行列を S を以下のように表す．

X̄n =パラメータ θ ∈ Θ を

n
n
1∑
1 ∑
Xj , S =(Xj − X̄n )(Xj − X̄n )⊤n j=1
n − 1 j=1

   ⊤ −1 
µ Σ µ
θ1
−1 



µ
θ = θ 2 = 1⊤
dΣ
−1
⊤
θ3
1d Σ 1dと定義する．ここで，以下の仮定を考える．

仮定 2.1.

(A1) n, d → ∞ かつ d/n → ρ ∈ (0, 1) とする．以後，これを (n, d)-asymptotic と呼ぶ．
1

(2.1)

(2.2)

(A2) 各 n, d を固定した下での θ を θ (n,d) ∈ Θ とする．このとき，ある θ (∞) = (α1 , α2 , α3 ) ∈ Θ ⊂R3 が存在して，(n, d)-asymptotic 下で以下が成り立つ．
θ (n,d) → θ (∞)(2.3)

このような条件の下で，次の統計量 θ̃ を定義する．

   ⊤ −1θ̃1X̄ S X̄nθ̃ = θ̃2  = 1d ⊤ S −1 X̄n 
1d ⊤ S −1 1d
θ̃3

(2.4)

このとき，次の結果が得られる．
定理 2.2. 仮定 2.1 下で次が成立する．
a.s.

θ̃1 →

1
ρ
α1 +
,
1−ρ
1−ρ

a.s.

θ̃2 →

1
α2 ,
1−ρ

a.s.

θ̃3 →

1
α3
1−ρ

(2.5)

n−d
d
⊤
系 2.3. θ̂ = n−dn ( n θ̃1 + n , θ̃2 , θ̃3 ) とする．このとき，仮定 2.1 の下で，次が成立する．
a.s.

θ̂ → θ (∞)(2.6)

この θ の一致推定量 θ̂ を用いて，有効フロンティアの一致推定量が得られる．
定理 2.4. 仮定 2.1 の下で

a.s.

sup ||ℓeﬀ (γ; θ̂) − ℓeﬀ (γ; θ (∞) )|| → 0

(2.7)

γ∈A

が成り立つ．
定理 2.5. 推定量 θ̂ = (θ̂1 , θ̂2 , θ̂3 )⊤ に対して仮定 2.1 の下で次が成立する．

√
Dn(θ̂ − θ) → N3 (0, Ω)
ここで，

(2.8)

 2
∗
2α + 4α1 + 2ρ1  1
2α1 α2
α22 + α1 α3 + α3Ω=1−ρ2α1 α32α2 α3


∗
∗ 
2α32

(2.9)

である．
定理 2.6. ℓeﬀ (γ; θ) = (σopt (γ; θ), µopt (γ; θ)) とし，ℓeﬀ は θ ∈ R3 まわりで微分可能であるとする．
このとき，

)
√ (
D
n ℓeﬀ (· ; θ̂) − ℓeﬀ (· ; θ) → N2 (0, ℓ̇eﬀ (· ; α)Ωℓ̇eﬀ (· ; α)⊤ )(2.10)

が成立する．ここで，(ℓ̇eﬀ (· ; θ) =

∂ℓeﬀ
∂θ

)⊤


=

である．

2

γ
γ2

−2γθ2 +1
θ3
−2γ 2 θ2
θ3

γθ22 −θ2θ32
γ 2 θ22 −1
2
θ3(2.11)

Rank tests for an ARMA model against other tv-ARMA models
Shunsuke Sakai and *Junichi Hirukawa
Abstract For a class of locally stationary process introduced by Dahlhaus, we apply the idea of the problem of
testing ARMA model against other non-stationary ARMA model. When testing the problem, we use linear serial rank
statistics and contiguity of LeCam’s notion. And then, if null hypothesis is white noise, under null and alternative, theasymptotic normality of the proposed statistics is established by using the locally asymptotic normal (LAN) property.
We incorporate the locally stationary phenomena in the testing problem.
Nonparametric methods have been developed for analysis of univariate and multivariate observations without the
distributional assumption (e.g., Gaussian assumption). The nonparametric procedures are even desired in the area oftime series analysis. Hallin et al. (1985) considered the systematic time series oriented study of testing for randomness
against ARMA alternatives. They proposed the statistic of the form
ST =
(T )

T
∑
1
(T )
(T )
(T )
cT (Rt , Rt−1 , . . . , Rt−p ),T − p t=p+1
(T )

(T )

(T )

where Rt denotes the rank of the observation Xt in an observed series X (T ) = (X1 , . . . , XT ) of length T , andcT (· · · ) is some given score function. These statistics are so-called linear serial rank statistics (of order p).

1

Locally asymptotic normality of time varying ARMA models

We consider the sequence of locally stationary ARMA(p1 , p2 ) (tv-ARMA) models(T )Xt,T +

p1
∑
k=1

∑
t
t
t
(T )
aθ(p2 +k) ( )Xt−k,T = bθ(0) ( )εt +bθ(k) ( )εt−k , t ∈ Z,
T
T
T
p2

k=1

{
}
where εt ; t ∈ Z is a discrete-time stationary white noise, i.e., a family of i.i.d. random variables with mean E[εt ] = 0,variance E[ε2t ] = 1, t ∈ Z and density function f . Let bθ(k) (u), aθ(p2 +k) (u) bebθ(k) (u) = θ(k) + (θ(k) − θk )(b̃θk (u) − 1), k = 0, . . . , p2 ,aθ(p2 +k) (u) = θ(p2 +k) + (θ(p2 +k) − θp2 +k )(ãθp2 +k (u) − 1), k = 1, . . . , p1 ,where the functions b̃θk (u), ãθp2 +k (u) are suﬃciently smooth bounded functions.
Consider the null hypothesis H(f ; θ) and the alternative hypothesis K(f ; θT ) as follows;hk(k)H(f ; θ) : θ(k) = θk (k = 0, . . . , p1 + p2 ), K(f ; θT ) : θ(k) = θT = θk + √ (k = 0, . . . , p1 + p2 ).
T′′

where θ = (θ0 , . . . , θp1 +p2 ) ∈ Θ ⊂ Rp1 +p2 +1 , h = (h0 , . . . , hp1 +p2 ) ∈ H ⊂ Rp1 +p2 +1 . Under the two hypotheses, thecoeﬃcient functions bθ(k) (u) and aθ(p2 +k) (u) are given by{{θk(H(f ; θ))θk(H(f ; θ))bθ(k) (u) =aθ(p2 +k) (u) =hhkkθk + √T b̃θk (u) (K(f ; θT )),θk + √T ãθp2 +k (u) (K(f ; θT )).
′Under null hypotheses H(f ; θ), for the Fisher information matrix Γ(θ) of Lemma 5, − 12 h Γ(θ)h can be written as
follows;∫p)
I(f ) 1 ∑(1 ′− h Γ(θ)h = −Ak (u) + Bk (u) + Ck (u) + Dk (u) du.
2
2
0
k=0

For the deﬁnitions of Ak (u), Bk (u), Ck (u) and Dk (u), refer to the master thesis.
(T )In the special case of the white noise null hypothesis H0 : Xt,T = θ0 εt , we can obtain)∫ ( pI(f ) 1 ∑ 2
1 ′dk (u) du,− h Γ(θ)h = −
2
2
0
k=0

where
√
(
)2

b̃θ0 (u)h01
(E[ε2t ϕ2 (εt )] − 1)

I(f
)
θ

0


dk (u) = ãθp +k (u)hp2 +k − θ1 b̃θk (u)hk
2
0



ãθp2 +k (u)hp2 +k 1− θ0 b̃θk (u)hk

k=0
1 ≤ k ≤ min(p1 , p2 )p2 < k ≤ p1 if p2 < p1
p1 < k ≤ p2 if p1 < p2 .

Under null hypothesis H(f ; θ), we derive the central sequence ∆T (θ) of Lemma 5 for locally stationary ARMA
model
T {
t−1
}
∑
∑
(t)
(t)
1
h′ ∆T (θ) = √ϕ(εt )
Dk′
Xt−k,T − b̃θ0
(1 + ϕ(εt )εt )h0 + op (1).
TT
T θ0 t=p+1
k=1( )For the deﬁnition of Dk′ Tt , refer to the master thesis.
(T )In the special case of the white noise null hypothesis H0 : Xt,T = θ0 εt , we have
p
T {
}
∑
∑
(t)
(t)
′
1
h ∆T (θ) = √ϕ(εt )dk
Xt−k,T − b̃θ0
(1 + ϕ(εt )εt )h0 + op (1),TT
T θ0 t=p+1
k=1

where

2


(t)
(t)
ã
h
− 1 b̃h( t )  θp2 +k ( Tt ) p2 +k θ0 θk T k= ãθp2 +k T hp2 +k
dk

( )
T
 1
− θ0 b̃θk Tt hk

1 ≤ k ≤ min(p1 , p2 )p2 < k ≤ p1 if p2 < p1
p1 < k ≤ p2 if p1 < p2 .

Main results

Now, we attempt to derive the asymptotic distribution
linear)serial rank statistics S (T ) under K(f ; θT ). If we( 1/2 (Tof) the (TT (S− m ))(T )can derive the joint asymptotic normality ofunder H0 , then it will lead to the asymptotic
ΛT (θ, θT )normality of linear serial rank statistics S (T ) under K(f ; θT ), from LeCam’s third lemma (see, e.g. Hajek et al.
(T )(1999)). In what follows, we assume that the null hypothesis is white noise hypothesis H0 : Xt,T = θ0 εt .
1/2
From Hallin and Puri (1985), we can establish the asymptotic equivalence of (T − p) (S (T ) − m(T ) ) with S (T ) −(T )M . (For the deﬁnitions S (T ) and M (T ) , refer to the master thesis). Furthermore, we can show that the U-statistics
(T )
(T )
US − UM which are asymptotically equivalent to T −1/2 (S (T ) − M (T ) ).
′(T )We deﬁne the U-statistics UL associated with the central sequence L (T ) = h ∆T (θ) as follows. Then, we showthat this U-statistics is asymptotically equivalent to the central sequence T −1/2 L (T ) .
Consider the function Gt,T (y) (of (p + 1) arguments y = (y1 , . . . , yp+1 ) ∈ [0, 1]p+1 ) where
Gt,T (y) = ϕ(F −1 (y1 ))

p
∑
i=1

di

( t ) −1}b̃θ ( t ) {F (yi+1 ) − 0 T 1 + ϕ(F −1 (y1 ))F −1 (y1 ) h0 .
Tθ0Additionally, put
(T )
(T )
ΦL
(t1 ,...,tp+1 ),T (Yt1 , . . . , Ytp+1 ) =

p+1
1 ∑
(T )
Gt ,T (Ytj ),p + 1 j=1 jwhere Yt
= (Yt,1 , . . . , Yt,p+1 )′ = (F (Xt,T /θ0 ), . . . , F (Xt−p,T /θ0 ))′ , p + 1 ≤ t ≤ T . Then, the corresponding
U-statistic is
(
)−1
∑
T −p
(T )
(T )
(T )
−1/2
ΦL
L (T ) + op (T −1/2 ).
UL =(t1 ,...,tp+1 ),T (Yt1 , . . . , Ytp+1 ) = T
p+1
(T )

(T )

(T )

p+1≤t1 <···<tp+1 ≤T

From the above, we can use the theory of asymptotic normality about U-statistics.
Our purposethat the)following joint asymptotic normality;( is1/2showing
)
( 2
(
)
V
Q
0
T (S (T ) − m(T ) )(T )′′Under H0 ,covarianceis joint asymptotically normal, with mean,ΛT (θ, θT )Q h Γ(θ)h− 12 h Γ(θ)hwhere V 2 is given by Hallin et al. (1985), Proposition 4.2.
If we can derive this joint asymptotic normality, from LeCam’s third lemma, we can derive the following desire
results;
Under K(f ; θT ), T 1/2 (S (T ) − m(T ) ) is asymptotically normal, with mean Q and variance V 2 .
Therefore, we need to prove this joint asymptotic normality and so that calculate the covariance Q.

References[1] Hallin, M., Ingenbleek, J-Fr., and Puri, M. L. (1985). Linear Serial Rank Tests for Randomness Against Arma
Alternatives. Ann Statist. 13 1156-1181.
[2] Hirukawa, J. and Taniguchi, M. (2006). LAN theorem for non-Gaussian locally stationary processes and its applications. Journal of Statistical Planning and Inference 136 640-688.

疎な分割表上の分布星野伸明

2015 年 11 月 17 日

1

高次元小標本（疎）な分割表

公的統計や社会調査では、個体属性の調査結果が離散変数とみなせる。このような状況は、変数が p 個として、標本の n 個体が p 次元の分割表上で分布していると考えられる。公的統計の場
合、典型的な p のオーダーは 102 程度で、n のオーダーは十万 (106 ) なら大きい方である。この場
合は n ≫ p なので、一見、高次元小標本とは言えない。しかし変数がほぼカテゴリカルなので、実質的なデータの次元は p ではない。一つのカテゴリカル変数が c 分類されていれば、実質的に
c − 1 個の変数が存在すると考えられる。故にカテゴリカルな分割表の実質的な次元 d のオーダー
は、cp である。典型的な c のオーダーは 10 なので、d ≍ 10100 ≫ 106 ≍ n となり高次元小標本と言えよう。
通常の漸近論は、中心極限定理を使いたいので n → ∞ とする。しかし統計的開示制限では n が所与での度数分布を用いるため、そのような極限操作は不都合である。一方、個体属性を詳細に表す場合に興味があり、それは c を増やすことに他ならない。その場合実質的な次元 d は増加する。従って本報告では n を固定して d → ∞ という漸近論を考える。中心極限定理の代わりに小数法則（ポアソン分布:Po）を用いれば、極限分布の取り扱いは解析的に容易となる。

2

疎な分割表の漸近論

分割表の実質的な次元 d のオーダーはセル総数 J と等しいので、以下では J → ∞ とする。ま∑た第 j セルの度数を Fj,J 、寸法指標を Si,J := Jj=1 I(Fj,J = i), i ∈ N, で表す。そのベクトルは
S J := (S1,J , S2,J , . . .) である。総度数は NJ で表す。
本報告の漸近論は次の命題 1 を用いる。なお複合ポアソン分布 CP(µ, q) を確率母関数 G(z) =∑iexp(λ(g(z) − 1)), λ > 0, で定義する。ただし g(z) = ∞i=1 qi z は正の整数上の分布の確率母関数
である（qi を i の確率と見る）。
命題 1 確率変数 Fj,J , j ∈ [J], が独立に非負整数上の分布に従うとする。この時 J → ∞ の極限でd∞

(S J , NJ ) → ( ⊥⊥ Po(µqi ), CP(µ, q))
i=1

(1)

が成立することは以下の条件 (2) と (3) が成立することと同値である。

lim max P(Fj,J = i) = 0,

J→∞

j

1

i ∈ N,(2)

lim E(Si,J ) = µqi ,J→∞i ∈ N.

(3)

命題 1 は Hoshino (2012) の結果を書き直した。寸法指標 Si,J の周辺分布はポアソンの二項分布をする。条件 (2) と (3) は、ポアソンの二項分布で各成功確率は小さいが成功回数の期待値が（正値に）収束することを意味している。つまり極限 (1) で寸法指標がポアソン分布に従う理屈は、小数法則に他ならない。
命題 1 において NJ の極限分布が CP(µ, q) となることは、以下のように説明出来る。確率変数
ベクトル∞

S := (S1 , S2 , . . .) ∼ ⊥⊥ Po(µqi )i=1としよう。この時 Si の確率母関数は exp(µqi (z −1)) であり、iSi の確率母関数は exp(µqi (z i −1)) と
∑
∏
∑ i
i
なる。各 Si は独立なので N := ∞i=1 iSi の確率母関数はi exp(µqi (z −1)) = exp(µ( i z qi −1))∑となる。これは CP(µ, q) の確率母関数であり、NJ = i iSi,J の極限分布を与えている。この独立なポアソン分布に従う確率変数の整数倍の和としての複合ポアソン分布の表現は “composed
Poisson” と呼ばれる
ここで S の N 所与での条件付き分布を考えよう。全ベル多項式を用いれば xi = i! qi の時P(N = n) =exp(−µ)Bn (µx1 , . . . , µxn ),n!

n ∈ N0 ,と書ける。N ∼ CP(µ, q) だったので、複合ポアソン分布の確率関数はベル多項式を用いて表現可
能ということになる。結局∏
n!µk ni=1 qi si s1i !
P(S = s|N = n) =,(4)
Bn (µx1 , . . . , µxn )∑ただし xi = i! qi , k = ni=1 si ということになる。
分布 (4) は CCP 分布の J → ∞ における極限分布である。従って分布 (4) を母数 (µ, q) の Limiting
Conditional Compound Poisson (LCCP) 分布と呼ぶ。値域を見れば、LCCP 分布は自然数の順序の着かない確率分割族であることが分かる。
LCCP 分布族は Ewens 分布など解析的に取り扱いが容易な例を含み、族としての性質も議論できる。これまでに報告者は周辺分布やモメントを明らかにしている。

参考文献[1] Hoshino, N. (2012). Random partitioning over a sparse contingency table. Annals of theInstitute of Statistical Mathematics, 64, 457–474.

2

Some multivariate asymmetric kernel estimators for density functions
柿沢 佳秀 (北大経済)1. はじめに カーネル密度推定量の漸近的性質はよく知られている (Silverman (1986),
Wand and Jones (1995)) が, 推定される密度 f の台が Rd でない (典型的に [0, 1]d や[0, ∞)d の台) のときに, 境界バイアス問題がある. カーネル密度推定の漸近論は, その




1
x −v
背後にある積分近似 Rd hd K h f (v) dv = Rd K(u)f (x − hu) du ≈ f (x) (h → 0)
が本質的 (ただし, K(−v) = K(v), Rd K(v) dv = 1 を仮定) で, 有界閉領域や正象限の境界付近ではこの近似が成り立たない. このような境界バイアス問題への対処として,境界カーネル法が主に１次元データに対し議論されてきた (Jones (1993)). 一方, 推定
密度 f が区間 [0, 1]/半無限区間 [0, ∞) の台を持つ場合に, Chen (1999,2000) がベータ/
ガンマカーネル密度推定量を提案して以来, 最近１０年間において『非対称カーネル法(AK 推定法)』が脚光を浴びてきた. 以下, 非負データを扱う.
2. １次元の AK 推定 １次元の非負データ X1 , . . . , Xn ∼ iid f に基づいた (点 x ∈
[0, ∞) = R+ での) AK 密度推定量は, n−1 ni=1 kx,b(Xi ) と定義される1 . ここで, kx,b(·)
は f と同じ台 R+ を持つカーネルとし, b = b(n) > 0 は n → ∞ のとき適当な速さ
で 0 に収束する平滑化パラメータである. このような先行研究として, ガンマカーネル(Chen (2000)), Birnbaum–Saunders (BS)/対数正規カーネル (Jin and Kawczak (2003)),
逆ガウス/相反逆ガウスカーネル (Scaillet (2004)), 逆ガンマカーネル (Koul and Song(2013)), 一般化 BS カーネル族 (Marchant et al. (2013)), 歪 BS カーネル族 (Saulo etal. (2013)), 一般化逆ガウスカーネル族 (Igarashi and Kakizawa (2014)), 一般化ガン
マカーネル族 (Hirukawa and Sakudo (2015)), 重み付き対数正規カーネル族 (Igarashi(2015; to appear)) がある.
3. 多次元の AK 推定 d 次元の非負データ X1 , . . . , Xn ∼ iid f (ただし, Xi =(Xi1 , . . . , Xid ) とする) に基づいて, f (x), x = (x1 , . . . , xd ) ∈ Rd+ を推定したいとする.
3.1. (積型の)AK 推定 カーネル kx,b(·) を 2 節のように選び, その積として d 次元の
カーネルを作れば, このときの (積型の)AK 密度推定量は n−1 ni=1 dj=1 kxj ,bj (Xij ) と定義される. ここで, bj = bj (n) ∝ b(n) とする (Bouezmarni and Rombouts (2010) が,Chen (1999,2000) のベータ/ガンマカーネルによる (積型の)AK 密度推定量を考察).
3.2. (非積型の)AK 推定 Koul and Song (2013) に述べられた注意点も一理があり,
多変量ガンマ分布 (Kotz et al. (2000) のテキストからその定義自体も複数存在, かつ，無限級数を伴う定義になっている) を直接, 密度推定に応用するのは無理があろうことは致し方ない (勿論, 理論的な関心は残されている).
一方, 要素毎に適当な変換を施して多変量分布が構成される場合, 計算機の面からも解析的にも扱いが容易になることが期待できる. 典型例として, (重み付き) 多変量対数
正規分布があり, これに基づいた (非積型の)AK 密度推定量を五十嵐 (2015; 統計関連
1

この形式は “重み付き推定量” (例えば, Siliverman (1986, Subsection 2.9: General weight function
estimators)) であるが, Chen (1999,2000) 以降は “非対称カーネル推定量”がキーワードである.

学会連合大会) が報告している. 本報告では, 多変量 BS 分布を (非積型の) カーネルとして用いることに焦点をおき, その拡張についても議論した (なお, 対数正規分布と BS
分布との関係も触れた).

4. 多変量 BS カーネル密度推定量とその拡張
4.1. 準備 (分布論) １次元の確率変数 ξ が密度関数k(s; α, β) = φ a(s/β)  A(s/β)ααβ(s ∈ R+ )(BS)(Birnbaum and Saunders (1969)) を持つとき, ξ は BS 分布 BS(α, β) に従うといい
,


1
1
1
1
1/2
ξ ∼ BS(α, β) と書く. ここに, α, β > 0 であり, a(t) = t − t1/2 , A(t) = 2 t1/2 + t3/2 ,2φ(u) = √12π exp − u2 . Kundu et al. (2010) は自然な拡張として２次元 BS(BBS) 分布BBS(α1 , α2 , β1 , β2 , ρ) を定義した (この d 次元化も同様に定義できるが, ここでは省略):
k(s; α1 , α2 , β1 , β2 , ρ) = φρ a(s /β )
1

α1

1

,

2
A(sj /βj )a(s2 /β2 )  
α2
αj βj
j=1



(s = (s1 , s2 ) ∈ R2+ ) .
(u2 +u2 −2ρu u )

1 2
1
1
2
ここに, α1 , α2 , β1 , β2 > 0, |ρ| < 1, φρ (u1 , u2 ) = 2π(1−ρ.
2 )1/2 exp −2(1−ρ2 )4.2. 密度推定への応用 想定したカーネル (ここでは多変量 BS, 及びその拡張) に付随するパラメータを推定位置に依存させることから, ノンパラメトリックな (可変的な)AK
密度推定量を構築できる. すなわち, 上記のパラメータ (α, β) に対し,b 1/2, βb,c (x) = x + bc (ただし, c > 0)x + bcを採用し, (非積型の) ２次元 BS カーネル密度推定量 (d 次元化も同様)
αb,c (x) =(BBS )

ρf(b1 ,b2 ),(c(x) =1 ,c2 )

1 n (BBSρ )k(Xi ; x) ,n i=1 (b1 ,b2 ),(c1 ,c2)

x = (x1 , x2 ) ∈ R2+を提案した (ρ = 0 のとき (積型の) 密度推定量に他ならない). さらに，漸近バイアス/
分散/平均積分２乗誤差/漸近正規性/強一致性の結果, 及び数値例を報告した.

4.3. 拡張 d 次元でも同様の議論があるので, ここでは１次元に限定する.
(I) BS 密度は, 標準正規密度 φ の変換分布である. その φ を他の対称密度にするとか,歪みを加えることもできる (BS 分布よりも柔軟なモデル化を扱った文献は膨大).
(II) 分布生成としての変換法の観点からは, おそらく BS 密度よりも対数正規密度11
exp − 2 (log s − μ)2
fLN (s) = √(s ∈ R+ )
(LN)
2σ
2πσs
がよく知られている. 変換による特徴付け『Z = a(ξ (BS) /β)/α ∼ N(0, 1)』の関数 a を⎧⎪⎨ 1 (tq − t−q ) , q = 0 ,[q]
2qa (t) =⎪⎩ log t ,q=0

(Yang (2006)) へと拡張することで両者を含む分布族になる.
これらの方向からの密度推定の研究は現在進行中である (進展があれば, 別の機会に報告したい).

事前密度のベイズ尤度とその含意

柳 本 武 美 統計数理研究所

 問題の設定
確率モデルを仮定して得られたデータ を評価する頻度論は、今日の統計学に確固たる基盤をなしている。一方、標本密度    を母数  の関数と見て尤度と見なすことがある。尤度に基づいた 尤度推論が優れていることが広く受け入れられている。対応するベイズモデルにおける尤度は事
前密度の確からしさである、との仮定の下に尤度の定義を試みている。尤度が導入できると様々な適用が期待できる。実際の適用では標本密度より事前密度の評価を行う場合が殆どだから、この報
告では確率密度の評価には言及しない。

 ベイズ尤度の定義
ベイズ尤度を予測子応じて定義する。その意味で複数のベイズ尤度を定義する。標本密度を   
事前密度を    とする。また、 を と同じ標本空間を上の将来 未観測 の標本として予測子は
一般に   と書く。

  予測子 により与えられ、     
ベイズモデルにおける予測子    は 
で定義される。ここで  は と同じ標本空間の変数であり、未来の観測値 あるいは未知観測値
とする。事前密度ではなく事後密度の下での期待値から導かれている。この予測子は、密度関数の空間の微分幾何学的接続の概念を用いて一般化できることを       
 が示した。彼らは一般的な接続係数  に対して様々な予測子を議論した。その両端である 混合と 混合が多くの場合に重要である。          後者の予測子が実際的あることを論じている。特に !"#$" %  との関連が深くなる。ベイズ尤度の関連では予測子を定義するためには、事後密度だけではなくて事前密度による混合も用いられる。結局以下のように四つの予測が議論の対象とする。

  
&

四つの予測子は次のように定義される

  ' (%)    * ½ +   '   &

,



-

 



 ' (%)     '   

ここで ¾ ' (% .  



*

¾ 







    ' (%

  &

 四つの尤度
予測子   からベイズ尤度を導くために     の不変な量を と    を用いて構成する。考え方は情報量規準の導出と同じである。その上で、改めて指数変換を施してベイズ尤度を導出する。この指針は多様なベイズ尤度の定義を導いてしまう意味で一意には定められない。

   & ./ 

で与えた予測しに対応して次の四つのベイズ尤度が定義される&

   '    +

 $ 



 $



 '

   +



 '

 



 '

     (%. , $ - $  (%      &  定義に見る尤度の特徴
尤度 $    は通常は周辺尤度と呼ばれる。問題は対応する予測子が事後密度ではなくて事前密度を用いて導出されている点である。尤度 $   を対数変換して  を乗じると 0%
ら   が導入した .1 の特定化版になる。特定化では予測子を前面に出して % 予測子を排除した点にあるからごく真っ当である。その結果、調整項  が常に正になる。.1 は 230 に搭載されるなど広く用いられているが、様々な議論がある 例えば 0% ら+  -。
多くの研究者は予測子      に基づいた尤度の導入を期待しているようである 例えば 
ら&  - 0% ら+  -。この予測子は対数変換すると不偏な量を見いだすのは容易ではない。その理由は、標本分布族にしばしば見られる積構造と 混合と折り合いの悪さに起因する。
これ以外のベイズ尤度は、小修正を除けばほぼ確定できる思われる。

 含意
この節では、ベイズ尤度が定義できたときの含意を議論する。

  経験ベイズモデルにおける超母数の推定尤度が定義できると最初の適用は経験ベイスモデルへの適用である。事前密度の評価は、超母数の評価と読み換えられるから、ごく自然に最大ベイズ尤度推定量により超母数を推定することにな
る。周辺尤度最大化と最小 .1 の比較についての研究は、二つのベイズ尤度を用いたときの性能の比較と見なされる。ベイズ尤度として一つの枠組みとして理解することにより高い視点からの比較が可能である。

 情報量規準
0% ら + % 45 による .1 は 1 のベイズ版として調整項 ペナルティ項 を変更したものである。既に述べたように、   特定版 .1 は   $   に一致する。そして、特に無情報事前密度が仮定された下では .1 が 1 と一致することがあることを指摘されている。
情報量規準 21 はサイズ１の標本相当の情報量を含む事前情報を仮定したときの周辺尤度で近似されることが知られている 。この事実は次のように見ることができる。情報量規準はベイズ尤度の特殊な場合である。

   はベイズ尤度が低い二つのベイズモデルの比較から生じている。より柔軟な二つのベイズモデルと繋ぐモデルの族を考える。その分布族に関してベイズ尤度を評価することにより現実的な推論が構成される。 推定量が意味あるものだとすれば、尤度比検定統計量も役立つと期待される。
$"#6 % (ランダム係数を持つ GMANOVA モデルの変数選択規準
広島大

理

若木

宏文

(l)

Yij (l = 1, . . . , g; i = 1, . . . , nl ; j = 1, . . . , p) を第 l 群に属する個体 i の tj 時点での観測値とし, 次の成長曲線モデルを考える.
(l)

(l)

Yij = bi +

q
∑

(l)

(l)

βk xk (tj ) + εij

(1)

k=1
(l)

(l)

ここで, xk (t) (k = 1, . . . , q) は基底関数, bi は個体変動を表わすランダム係数, εij は観
測誤差で,i.i.d.

(l)(l)

bi (i = 1, . . . , nl ; l = 1, . . . , g) ∼ N (β0 , τ 2 ),i.i.d.

(l)εij (l = 1, . . . , g; i = 1, . . . , nl ; j = 1, . . . , p) ∼ N (0, σ 2 )とする.
θ をすべての未知母数を並べたベクトル, Y をすべての観測変量を並べた確率行列とし, θ̂(Y) を Y による最尤推定量とすると, カルバック-ライブラーダイヴァージェンスに基づくモデルの良さは, Z を Y と独立に同じ分布に従う確率行列としてRisk(f ) = −2EY EZ [log f (Z; θ̂(Y))](2)

で与えられる. ここで, f (Y, θ) は, Y の同時確率密度関数である.

 (l)  (l) 
Yi1
εi1
x1 (t1 ) · · · ..  (l)  ..  ..
(l)Yi =  .  , εi =  .  , X =  .
(l)

Yip

 (1)
β1 · · ·
 ..
B= .
(1)

βq

···

(l)

(g) 

εip

x1 (tp ) · · ·xq (t1 ).. . 
xq (tp )

(1)
(g)
β1
bi = (bi , . . . , bi ),..  ,. (g)el = (0, . . . , 0, 1, 0, . . . , 0)TβqlT
とおき, 説明変数行列について rank X = q, 1T
p X = 0p と仮定する.

λ2 = σ 2 + pτ 2 ,

B = (β (1) , . . . , β (l) , . . . , β (g) )(l)

とおくと, σ 2 , λ2 , β0 , β (l) (l = 1, . . . , g) の最尤推定量は

1 T (l)1(l)βb0 =1p Y 1nl , βb(l) = (XT X)−1 XT Y(l) 1nl ,
nl p
n
 (
)l
S2
S1


,
if S1 ≤ (p − 1)S2n(p − 1) n2 b2(bσ ,λ ) = ()

S1 + S2 S1 + S2
,
if S1 > (p − 1)S2

npnp

で与えられる. ここで,g {}∑
1 T (l) T
(l) T
(l)
(l)
tr[(Ip − Gp )Y (Inl − Gnl )(Y ) ] + 1nl (Y ) (Ip − Gp − ΠX )Y 1nl
S1 =
nl
l=1

1 ∑ T (l)1p Y (Inl − Gnl )(Y(l) )T 1p ,
p l=1
g

S2 =

(l)

Y(l) = [Y1 , . . . , Yn(l)], Gnl =
l

1
1
T
T
−1 T
1 nl 1 T
nl , Gp = 1p 1p , ΠX = X(X X) X
nlpである. (2) を (−2)× 最大対数尤度 −2 log f (Y; θ̂(Y)) によって推定するときのバイア
スは[]S2 − n 2 λ 2
S1 − n 1 σ 2
ES1 ,S2+S1 {n(p − 1)}−1S2 n−1[{}]n1 σ 2n2 λ2n 1 σ 2 + n 2 λ2
+ ES1 ,S2 1{S1 >(p−1)S2 }+−,S1 {n(p − 1)}−1 S2 n−1 (S1 + S2 )(np)−1

n1 = n(p − 1) + gq, n2 = n − gと表される. ラプラス展開を用いてバイアスを評価することにより次の結果を得た.
定理(0, 1) 上の C 3 -級関数 H(x) の x = p0 = (p − 1)/p における微係数がH(p0 ) = 1, H ′ (p0 ) =2p2
gp(p − q − 1)
, H ′′ (p0 ) = −p−1p−1であるとする.

(
)
(
)
√
√
d1 = 2 −p + 2 − p2 − 2p + 2 , d2 = 2 −p + 2 + p2 − 2p + 2 ,c=2(p2 − p + 1) + g(3p2 − (7 + q)p + 4(q + 1))
p(p − 1)

とし,

[
]
(qg
+
1)(qg
+
2)
(g
+
1)(g
+
2)
b
AIC2 = −2 log f (Y; θ(Y))+ 2 {g(q + 1) + 2} ++n(p − 1) − qg − 2
n−g−2
{ (
}
)
}
{
np
S1
c
−
+ 1{ S1 >p−1+ d2 } H
+
,
1 S1
d1
S2
n
np − (q + 1)g − 2 { S2 >p−1+ n }S1 + S2nと定義するとE[AIC2 ] − Risk(f ) = O(n−3/2 )
を満たす.

[
]
AIC2 = −2 log f (Y; θ̂(Y)) + 2 {g(q + 1) + 2} − 1{S1 >(p−1)S2 } + Op (n−1/2 )であり, σ 2 , λ2 の最尤推定値がパラメータ空間の境界 {(σ 2 , λ2 ); σ 2 = λ2 > 0} 上にあるときは, 通常の AIC の罰則項よりパラメータ数をほぼ一つ分 (罰則項としては 2) 減らすべきことがわかった. また, 数値実験により O(n−1 ) までバイアスを補正した AIC2 の有効性が確認された.

Distribution Theory and Its Asymptotics in Shape Analysis(Joint Work with Professor Peter E. Jupp, University of St. Andrews, U. K.)
Yasuko Chikuse
Professor Emeritus, Kagawa University

1

INTRODUCTION

This paper is concerned with developing some families of distributions on the shapespaces. The “shape” of an object is the geometrical information that remains when allowance is made for changes in location, scale and orientation. The shape of a set X (m × k)
of k (not totally coincident) labelled points in the m-dimensional real space is transformedby the Helmert matrix (centred and scaled) to give a “pre-shape” Z (m × (k − 1)). The“(similarity) shape space” is then the coset space of the pre-shape space obtained by the
special orthogonal transformations U Z of Z for U in the special orthogonal group SO(m).
It is useful to identify distributions on the shape space with distributions of Z that areinvariant under the special orthogonal transformations U Z for U in SO(m). Here we canconstruct the “reﬂection shape space” obtained by the orthogonal transformations U Z for
U in the orthogonal group O(m), removing the eﬀect of both rotations and reﬂections. The
distributions on this reﬂection shape space can be identiﬁed with distributions of Z thatare invariant under the orthogonal transformations U Z for U in O(m).
The main parametric families of distributions which are concerned in this paper are
1. Shape Bingham distributions having densities proportional to exp[tr(AZ ′ Z)], and
2. Determinantal shape angular central Gaussian distributions having densities proportional to some power of the determinant of ZAZ ′ .
We consider the derivations of the normalizing constants of these distributions.
Also we can consider
3. the distributions without reﬂective symmetry, the volume Fisher-Bingham distributions, and the cardioid-type distributions.
Asymptotic behaviours of the distributions are considered for the large concentrations
and for the large labels.

REFERENCES[1] Chikuse, Yasuko and Jupp, Peter (2015). Some families of distributions on higher shapespaces, In “Geometry Driven Statistics”, Dryden and Kent, eds, Wiley, pp.206-217.
[2] Chikuse, Yasuko and Jupp, Peter (2004). A test of uniformity on shape spaces, J.
Multivariate Analysis 88, 163-176.
1[3] Chikuse, Yasuko (2003). “Statistics on Special Manoifolds”, Lecture Notes in Statistics,
Vol. 174, Springer.
[4] Mardia, K. and Jupp, P. (2000). “Directional Statistics”, Wiley, New York.
[5] Mardia, K. and Dryden, I. (1998). “Statistical Shape Analysis”, Wiley.
[6] Goodall, C. and Mardia, K. (1993). Multivariate aspects od shape theory, Annals of
Statistics, 21, 848-566.
[7] Davis, A. W, (1980). Invariant polynomials with two matrix arguments, extending thezonal polynomials, In “Multivariate Analysis, vol. V”, Krishnaiah, ed, pp. 287-299.
[8] Andrews, G., Askey, R. and Roy, R. (1999). “Special Functions”, Encyclopedia of Mathematics and Its Applications, Vol. 71, Cambridge University Press.
[9] Azzalini, A. (2013). “The Skew-Normal and Related Families”, IMS Monographs Series,
Cambridge University Press.

2修正コレスキー分解を用いた分散共分散行列の縮小推定について
広島大学 茶之原直人

中川智之

若木宏文

1 導入
本研究は共分散行列の推定に関する研究である. 従来の共分散行列の推定方法としては標本共分散行列が用いられ, 大標本下においては不偏性・正定値性が保証される面で良い推定が可能である. しかし, 次元数が標本数と比べて小さくない場合 (以下高次元と呼ぶ) においては, 最大(最小) 固有値のバイアスが増加する (Johnstone(2001)）ことや, 正定値性が保証されないことが問題視されている. したがって, 高次元下における共分散行列のより良い推定量が望まれており, 数多くの手法が提案されている (詳しくは Pourahmadi(2011) で多く紹介されている). 中でもPourahmadi(1999) による修正コレスキー分解を用いた回帰問題への帰着は, 推定量が必ず正定値性が成立すること, LASSO(Tibushirani(1996)) などの線形回帰モデルのテクニックが応用できることから, Huang et al(2006), Chang and Tsay(2010), Levina, Rothman and Zhu(2008) な
ど修正コレスキー分解を用いた罰則付き最尤推定法が提案されている. しかし多くの手法では共分散構造 (モデル) を限定した上での理論的性質しか分かっておらず, 構造が未知の場合の良さは保証されていない. そこで本研究では, 共分散構造を正定値性のみの仮定の下で, 修正コレスキー分解を用いた回帰係数の縮小型推定量の提案, 大標本下におけるリスクの漸近近似から理論的性質の導出を行った. 漸近近似の結果から, 罰則パラメータの選択法を提案するが, 罰則パラメータに真の共分散行列が含まれているため, その推定量について数値実験を行った.

2 提案手法
i.i.d.

仮定 1. Y = (y 1 , · · · , y n )′ , y i ∼ Np (0, Σ), Σ = (σij ) : 正定値
仮定 1 の下で, Σの推定を行いたいとする. 本研究では修正コレスキー分解 Pourahmadi(1999)を用いた推定量として, 回帰係数 ϕj に対する縮小推定を用いた推定量Σ̂P LSE を提案する.

()−1 ′ϕ̂j P LSE = Y ′j−1 Y j−1 + λj Ij−1
Y j−1 y (j)σ̂j2 P LSE = n1 (y (j) − Y j−1 ϕ̂j P LSE )′ (y (j) − Y j−1 ϕ̂j P LSE )j = 2, · · · , p(1)

ここで, y (k) = (y1,k , · · · , yn,k ), Y k = (y (1) , · · · , y (k) ), σ̂12 P LSE = y ′(1) y (1) /n, λj ≥ 0であり, λj = O(1) とする. 回帰残差平方和に Ridge 型罰則 (Hoerl and Kennard(1974)) を課したもので, 高次元下 (n < p) においても, λj > 0 であれば残差平方和は 0 にならないため, どのような状況下においても正定値な推定量が構成できる.

3 リスクの漸近近似
KL ロスに基づくリスクの漸近近似を用いて罰則パラメータの選択を考える. 罰則パラメータとして, 別々の罰則 λ = (λ2 , · · · , λp )′ と共通の罰則 (λ2 , · · · , λp )′ = (λ, · · · , λ)′ の 2 通りを (1) へ代入したときの KL リスクの最小化を行った.
系 1. n > p, n は十分大きい, p : 固定, (σ1k , · · · , σk−1,k )′ ̸= 0(∀ k = 2, · · · , p)(※) の下
で, O(n−2 ) までのリスクを最小とするパラメータ λ∗ , λ∗ とそのリスクは次である.

[
] p
{p (
}
)
5
E KL(Σ, Σ̂P LSE (λ∗ )) = (p + 1)n−1 +10p2 + 21p + 13 − f (Σ) n−2 + O(n− 2 )
2 {
12
}
p
∑ (
) σkk·k−1
σkk·k−1−1 2
f (Σ) =
trΣk−1, λk∗ = trΣ−1, k = 2, · · · , p
k−1 ′ −1
′ −1
ϕk Σk−1 ϕk
ϕk Σk−1 ϕk
k=2
{p (
}
[
] p
)
5
10p2 + 21p + 13 − g(Σ) + O(n− 2 )E KL(Σ, Σ̂P LSE (λ∗ )) = (p + 1)n−1 +
2
12
(∑p
)
∑p
2
−1
−1
k=2 trΣk−1
k=2 trΣk−1,
λ=g(Σ) =∗∑p ϕ′k Σ−1∑p ϕ′k Σ−1
k−1 ϕkk−1 ϕkk=2

k=2

σkk·k−1

σkk·k−1

ここで, Σj = (σkl ), 1 ≤ k, l ≤ j, σjj·j−1 = σjj − ϕ′j Σj−1 ϕj , j = 2, · · · , p である.

(※) 例えば (σ1k , · · · , σk−1,k )′ = 0 とする. このときλk∗ は, ϕk = 0 であるため, 分母が 0 となりλk∗ = ∞ となる. これはλk = O(1) であることを満たしていないので, このような共分散構
造下では漸近近似の精度は保証されない. しかしながら, λk∗ = ∞ である場合ϕ̂k P LSE = 0 となるため, 正しい推定 (罰則) を行っていると言える.

参考文献[1] Chang, C. and Tsay, R. S. (2010). Estimation of covariance matrix via the sparse
Cholesky factor with lasso. J. Statist. P lann. Inf erence 140 3858-3873.
[2] Hoerl, A. E. and Kennard, R. W. (1970). Ridge Regression. Applications to nonorthogonal problems. T echnometrics. 12.
[3] Huang, J. Z., Liu, N., Pourahmadi, M. and Liu, L. (2006). Covariance matrix selectionand estimation via penalised normal likelihood. Biometrika 93 85-98.
[4] Johnstone, I. M. (2001). On the distribution of the largest eigenvalue in principal components analysis. Ann. Statist. 29 295-327.
[5] Levina, E., Rothman, A. and Zhu, J. (2008). Sparse estimation of large covariance matrices via a nested Lasso penalty. Ann. Appl. Statist. 2 245-263.
[6] Pourahmadi, M. (2011). Covariance estimation: The GLM and Regularization perspectives. Statist. Sci. 26 369-387.
[7] Pourahmadi, M. (1999). Joint mean-covariance models with applications to longitudinal
data: Unconstrained parameterisation. Biometrika 86 677-690.
[8] Tibshirani,R.

(1996).

Regression

J. Roy. Statist. Soc. Ser. B 58 267-288.

shrinkageandselection

via

the

lasso.

順序カテゴリ正方分割表における対角指数対称モデルについて
渋谷
生亀
富澤

明
清貴
貞男

(東京理科大学大学院 理工学研究科)(東京理科大学 理工学部)(東京理科大学 理工学部)

1. はじめにR×R 順序カテゴリ正方分割表において, (i, j) セル確率を pij とする (i = 1, . . . , R; j =1, . . . , R). 対角指数対称 (DES) モデルは次のように定義される (Tomizawa, 1992):{
δ i+j d|j−i| (i ̸= j),
pij =
ψii
(i = j).
DES モデルは, セル確率が主対角線に対して対称であり, さらに主対角線に沿って指数的に変化する構造をもつ. 準対角指数対称 (QDES) モデルは次のように定義される (Iki,Yamamoto and Tomizawa, 2014):{
αi β j d|j−i| (i ̸= j),
pij =
ψii
(i = j).
特に α = β とおいたこのモデルは DES モデルである. 行変数を X, 列変数を Y とし, 平均一致 (ME) モデルを E(X) = E(Y ) として定義する. Iki et al. (2014) は定理「DES モ
デルが成り立つための必要十分条件は, QDES モデルと ME モデルの両方が成り立つことである」を与えた.
QDES モデルは, 潜在分布として周辺分散が等しい 2 変量正規分布が想定される場合
においてよく適合する (Iki et al., 2014). 本講演では, 潜在分布として周辺分散の等分散性を仮定しない 2 変量正規分布が想定される場合において, よく適合するモデルを提案する. さらに, DES モデルと QDES モデルの一般化を考える.

2. モデルの提案
モデルを次のように提案する:
{

pij =

2

2

δ1i+j δ2i +j d|j−i| (i ̸= j),ψii(i = j).

特に δ2 = 1 とおいたこのモデルは DES モデルである. よって, このモデルを拡張対角指数対称 (EDES) モデルと呼ぶことにする. さらに, モデルを次のように提案する:{
22
α1i α2i β1j β2j d|j−i| (i ̸= j),
pij =
ψii
(i = j).
特に α2 = β2 = 1 とおいたこのモデルは QDES モデルである. よって, このモデルを拡
張準対角指数対称 (EQDES) モデルと呼ぶことにする. また, α1 = β1 かつ α2 = β2 = 1
とおいたこのモデルは DES モデル, α1 = β1 かつ α2 = β2 とおいたこのモデルは EDES
モデルである.

3. モデルの分解
平均分散一致 (MVE) モデルを E(X) = E(Y ) かつ Var(X) = Var(Y ) として定義す
る. 次の定理を得る:
定理 1. EDES モデルが成り立つための必要十分条件は, EQDES モデルと MVE モデルの両方が成り立つことである.
モデルの適合度を検定するための尤度比カイ二乗統計量を G2 とする.
定理 2. EDES モデルに対する G2 は, EQDES モデルに対する G2 と MVE モデルに対する G2 との和に漸近的に同等である.

4. モデルの一般化
DES モデルと QDES モデルの一般化を考える. 任意に与えられた k (k = 1, . . . , R−1)に対して, モデルを次のように提案する:{ ∏kit +j td|j−i| (i ̸= j),t=1 δtpij =ψii(i = j).
k = 1 のときは DES モデルであり, k = 2 のときは EDES モデルである. すなわち, このモデルは DES モデルを一般化させたモデルである. このモデルを k 次対角指数対称(DES(k)) モデルと呼ぶことにする. さらに, 任意に与えられた k (k = 1, . . . , R − 1) に対して, モデルを次のように提案する:
{ ∏
t
k
it j
t=1 αt βt d|j−i| (i ̸= j),
pij =
ψii
(i = j).
k = 1 のときは QDES モデルであり, k = 2 のときは EQDES モデルである. すなわ
ち, このモデルは QDES モデルを一般化させたモデルである. このモデルを k 次準対角
指数対称 (QDES(k)) モデルと呼ぶことにする. また, αt = βt (t = 1, . . . , k) とおいた
QDES(k) モデルは DES(k) モデルである.
任意に与えられた k (k = 1, . . . , R − 1) に対して, k 次周辺積率一致 (MME(k)) モデ
ルを E(X t ) = E(Y t ) (t = 1, . . . , k) として定義する. 次の定理を得る:
定理 3. 任意に与えられた k (k = 1, . . . , R − 1) に対して, DES(k) モデルが成り立
つための必要十分条件は, QDES(k) モデルと MME(k) モデルの両方が成り立つことである.
定理 4. 任意に与えられた k (k = 1, . . . , R − 1) に対して, DES(k) モデルに対するG は, QDES(k) モデルに対する G2 と MME(k) モデルに対する G2 との和に漸近的に同等である.
2

参考文献[1] Iki, K., Yamamoto, K. and Tomizawa, S. (2014). Statistics and Probability Letters, 92,
33-38.
[2] Tomizawa, S. (1992). Journal of Applied Statistics, 19, 509-512.

パラレルプロファイルモデルにおけるランダムエフェクト共分散構造に対する尤度比
統計量の高次元漸近展開
広島大学大学院 理学研究科

稲津佑

(g)

パラレルプロファイルモデルは多変量解析の場においてよく用いられている重要なモデルである. xj

を第 g グループの j 番目の個体から得られた p 次元目的変数ベクトルとする. ただし, j = 1, . . . , Ng ,g = 1, . . . , k であり, N を N = N1 + · · · + Nk とする. また, N − p − k − 1 > 0 は常に成り立っている
(1)

(k)

ものとし, k は 2 以上の固定された自然数とする. ここで, x1 , . . . , xNk は互いに独立に以下の分布,(g)xj∼ Np (µ(g) , Σ), µ(g) = δg 1p + µ, µ = (µ1 , . . . , µp )′ , δk = 0,

(1)に従っていると仮定する. このとき, モデル (1) をパラレルプロファイルモデルと呼ぶ. このモデルに対して, Yokoyama & Fujikoshi (1993) は Σ に以下のランダムエフェクト共分散構造を仮定した.

Σ = λ2 1p 1′p + σ 2 Ip , λ2 ≥ 0, σ 2 > 0.

(2)

この Σ に対する強い仮定により, 彼らは平均パラメータ δg , µ に対するより有効な推定量, もしくはより強力な検定を導出することに成功している. とはいえ, これらの結果は Σ が (2) の構造をしているときに成り立つ結果であるから, Σ がランダムエフェクト共分散構造をしているかどうかを検定することが重要となる.

Yokoyama (1995) は (2) の妥当性を評価するための次の仮説
H0 : Σ = λ2 1p 1′p + σ 2 Ip , λ2 ≥ 0, σ 2 > 0 v.s. H1 : not H0 ,(3)について考察しており, (3) を検定するための尤度比検定統計量 −2 log Λ を導出している. Yokoyama(1995) は −2 log Λ の帰無分布の漸近分布を, Inatsu & Wakaki (2015) は −2 log Λ の帰無分布の漸近展
開公式をそれぞれ導出しており, これらを用いることで (3) を検定することが可能となる.
しかしながら, Yokoyama (1995) あるいは Inatsu & Wakaki (2015) の結果はいずれも, 次元数 p を固定した下で N を無限大に飛ばす大標本漸近理論に基づいている. そのため, 次元数 p も N と同じくらい大きい高次元の枠組みでは, 近似精度は極めて悪い. そこで, 本研究では以下に仮定する高次元の枠組み,(HD1) N → ∞ かつ p → ∞.
(HD2) (N − p) → ∞.の下, −2 log Λ を基準化した統計量 V の高次元漸近展開公式を導出した.
パラメータ λ2 > 0 の場合と λ2 = 0 の 2 通りの場合において, 高次元漸近展開公式を導出した. まず

λ2 > 0 の場合, V の帰無分布の高次元漸近展開は以下で与えられる.
{p
if p/N → ϱ ∈ [0, 1)−3P(V ≤ c) = Q2 (c) + O(ϑ ), ϑ =.
N − p if p/N → 1, N − p → ∞ただし, Q2 (c) は

( (4)
)}
{ (3)
κ̃
(κ̃(3) )2κ̃
h2 (c) +h3 (c) +h5 (c),
Q2 (c) = Φ(c) − ϕ(c)
6
24
72
1

で定義され, Φ(c), ϕ(c) はそれぞれ標準正規分布の分布関数と密度関数である. また, hr (x) は

h2 (x) = x2 − 1, h3 (x) = x3 − 3x, h5 (x) = x5 − 10x3 + 15x,である. また, κ̃(s) = κ(s) {κ(2) }−s/2 (s ≥ 3) であり,

{
(
)
(
)}
N − (p − 1) − kN −kκ(s) = (−1)s ψ (s−1)− ψ (s−1)
2
2
{
(
)
(
)}
p−2
∑
N −l−1
N −1
s
(s−1)
(s−1)
+ (−1)
ψ
−ψ
2
2
l=1
{
(
)
(
)}
N −1
(N − 1)(p − 1)s(s−1)
s (s−1)+ (−1) (p − 1)ψ
− (p − 1) ψ
2
2

(s = 1, 2, . . .),となる. ψ (s) (α) は

(
ψ

(s)

(α) =

d
dα

)s+1
log Γ(α) =


)
∞ (
∑

1
1

 −CE +−,


1+k k+α
k=0
∞
∑







k=0

(−1)s+1 s!
,(k + α)s+1(s = 0),,(s = 1, 2, . . .),で表され, CE はオイラー定数である.
一方, λ2 = 0 の場合, V の帰無分布の高次元漸近展開は以下で与えられる.

P(V ≤ c) = Q2 (c) +

2
(j)
∑
Q (c)
2

j=1

j!

(−L−1 )j Hj + O(ϑ∗ ),

ただし,

(
(j)

Q2 (c) =ddcϑ∗ = max{p−1 N −3/2 , p−3 , ϑ−3 }.

)jQ2 (c),

L=N√
κ(2) ,

であり, Hj は以下で与えられる.

(
Hj = 2jΓ(1/2 + j)
J1 Γ(2/2 + j) J2 Γ(3/2 + j)J0+√+
Γ(1/2)N Γ(3/2)
N Γ(2/2)).

係数 J0 , J1 , J2 はそれぞれ以下で与えられる.

1
J0 =
2
J2 =

{
}
3k 2 (p − 1)2 + 5p2 − 2p − 13k(p − 1) − p − 1√
1−, J1 =,
12N p(p − 1)6 πp(p − 1)

1 3k 2 (p − 1)2 + 5p2 − 2p − 1.
212p(p − 1)

References
Inatsu, Y. and Wakaki, H. (2015). Asymptotic expansions of the null distribution of the LR test
statistic for random-eﬀects covariance structure in a parallel profile model. TR 15–14, Statistical Research Group, Hiroshima University, Hiroshima.
Yokoyama, T. (1995). LR test for random-eﬀects covariance structure in a parallel profile model.
Ann. Inst. Statist. Math., 47, 309–320.
Yokoyama, T. and Fujikoshi, Y. (1993). A parallel profile model with random-eﬀects covariancestructure. J. Japan Statist. Soc., 23, 83–89.

2

マルコフ連鎖モンテカルロ法の高次元漸近論

大阪大学大学院基礎工学研究科，金融・保険教育研究センター，JST CREST

鎌谷研吾

マルコフ連鎖モンテカルロ (MCMC; Metropolis et al. [1953]) 法は終戦後のロスアラモス研究所で開発された積分の近似計算法の一つである．Pd (dx) = pd (x)dx は Rd 上の確率分布とする．最も基本的な MCMC 法であるランダムウォーク型メトロポリス (RWM) 法は次の手続きで定義された．ある x ∈ Rd を出発点として以下を繰り返す：

Random walk Metropolis algorithm• x∗ ← x + w, w ∼ Nd (0, σ 2 Id ).
• x ← x∗ with probability α(x, x∗ ), ただし，{}pd (x∗ )α(x, x∗ ) = min 1,.
pd (x)d
上記のステップにより d 次元の確率変数列 {Xm}m が定義される．このマルコフ連鎖は Pd (dx) を不変分布として持つように出来ていて，例えば pd (x) > 0 (x ∈ Rd ) であればエルゴード的になり，下記右辺の積分が存在する限り大数の法則

1 ∑
d
f (Xm
)→
M m=1
M

∫
f (x)Pd (dx) a.s. (M → ∞)(1)が成り立つ．したがって右辺の積分の近似値として左辺の有限和を用いることが正当化されるわけである．

MCMC 法は Gelfand and Smith [1990] を契機としてベイズ統計学で利用されるようになり，特に高次元での積分計算には欠かせなくなった．状態空間の次元の高さはしかし，MCMC 法のパフォーマンスを低下させるd
主要因でもある．高次元であればあるほど pd (x) の計算負荷が増加する上，{Xm}m の動きも制限されるからである．MCMC 法の高次元漸近論はしばしばこれらの次元の影響を見積もる際に有用な情報を提供する．

MCMC 法の高次元漸近論はもともと量子力学の分野で研究されていた (Creutz [1988])．しかしパラメータ次元 d の増加に対する採択率 α(x, x∗ ) の安定性の解析であり，MCMC 法のパフォーマンスとのつながりは明瞭とは言えなかった．転機は Roberts et al. [1997] らによるベイズ統計学的な解釈で，彼らの研究により，高次元漸近論が MCMC 法のチューニングパラメータの調整の Criterion として使われるようになる，実用面でのインパクトをもたらした．ここで彼らの結果を端的に述べる．
興味のある分布は次の形であるとする*1 ．

Pd (dx) =d∏f (xi )dxi (x = (x1 , . . . , xd )).

i=1
d
RWM 法により d 次元の系列 {Xm}m が定義され，各 d, m でdddXm= (Xm,1, . . . , Xm,m).

*1 この仮定は一般的のように見えるかもしれないが，ベイズ統計学での事後分布に対応する量であるから，現実から乖離したものであり，またおおよそ正規分布に近似されることも意味する．

dと書けるとする．ここで時間スケーリング：t 7→ [dt] (dt ∈ R の整数部分) を定義し，X[dt]とすることにd
より連続時間の確率過程として {Xm}m を捉える．また任意の集合 E ⊂ {1, . . . , d} に対し射影 πE (x) =d(xi )i∈E (x = (x1 , . . . , xd ) を定める．その上で Ytd := π{1} (X[dt]). と置く．提案分散のスケーリングとして，

σ 2 = l2 /d と定める．
Theorem 1 (Roberts et al. [1997])
正則条件かつ定常性のもと，Y d ⇒ Y wheredYt = a(Yt )dt + b(Yt )dWt ,where(log f )′ (y)h(l), b(y)2 = h(l)
2(
√ )
)2
∫ ( ′
l I
f (x)
2
h(l) =2l Φ −, I=
f (x)dx < ∞,2f (x)a(y) =

Roberts et al. [1997] らはこの定理に実用的解釈を与え，それが大きなインパクトを与えた．本発表と関連する点は，上述の定理において t 7→ [dt] とした事実である．これにより，RWM 法は次元が高くなるつれ．d のオーダーで繰り返し回数を増やしていかなければいけないことがわかる．
本発表では様々な MCMC 法の高次元漸近論を考え，手法に比較や条件の比較を行うことで，高次元における有用な MCMC 法のデザインを考えることが目的である．主に次のアルゴリズム，MpCN 法を解析する．

• r ∼ Gamma(d/2, ∥x − µ∥2 /2),• x∗ ← µ + ρ1/2 (x − µ) + (1 − ρ)1/2 r−1/2 w, w ∼ Nd (0, σ 2 Id ),• x ← x∗ with probability α(x, x∗ ), ただし，{p(y)∥x − µ∥−dα(x, y) = min 1,p(x)∥y − µ∥−d}.

本発表は Kamatani [2014a,b] に基づいている．また，本研究は JSPS 科研費 24740062 の助成を受けている．

参考文献
Michael Creutz. Global monte carlo algorithms for many-fermion systems. Phys. Rev. D, 38:1228–1238, Aug
1988. doi: 10.1103/PhysRevD.38.1228.
Alan E. Gelfand and Adrian F. M. Smith. Sampling-based approaches to calculating marginal densities. J. Amer.
Statist. Assoc., 85(410):398, 1990.
Kengo Kamatani. Rate optimality of Random walk Metropolis algorithm in high-dimension with heavy-tailed
target distribution. Arxiv, 2014a. URL http://arxiv.org/abs/1406.5392.
Kengo Kamatani. Efficient strategy for the Markov chain Monte Carlo in high-dimension with heavy-tailed target
probability distribution. Arxiv, 2014b. URL http://arxiv.org/abs/1412.6231.
N. Metropolis, W. A. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller. Equations of state calculations byfast computing machines. Journal of Chemical Physics, 21:1087–1092, 1953.
Gareth O. Roberts, Andrew Gelman, and Walter R. Gilks.
dom walk Metropolis algorithms.
10.1214/aoap/1034625254.

Weak convergence and optimal scaling of ran-Ann. Appl. Probab., 7(1):110–120, 1997.

ISSN 1050-5164.

doi: