科研費シンポジウム「多様な分野における統計科学の総合的研究」
日時：2017 年 11 月 17 日（金）〜11 月 19 日（日）(Date: Friday, 17, November – Sunday, 19, November, 2017)
場所：コープシティ花園 4F ガレッソホール a(Place: Coop City Hanazono 4F Garesso Hall a)（TEL: 025-248-7511）
科学研究費・基盤研究（A) （課題番号：15H01678）「大規模複雑データの理論と方法論の総合的研究」
研究代表者：青嶋 誠（筑波大学）
Makoto Aoshima (University of Tsukuba)
開催責任者：蛭川 潤一（新潟大学）
Junichi Hirukawa (Niigata Universiry)Program
Friday, 17, November
Reception 13:00-13:10
Opening 13:10-13:15 Junichi HIRUKAWA (Niigata University)
Afternoon Session I (in Japanese) 13:15-15:15
Chair Takuma Yoshida (Kagoshima University)
1. 13:15-13:55 阿部 文貴 (Abe Fumitaka), 白石 博 (Shiraishi Hiroshi)
慶應義塾大学理工学研究科 (Keio University Graduate School of Science and Technology)
慶應義塾大学理工学部数理科学科 (Keio University Department of Mathematics)「長期記憶性を持った高次元ポートフォリオの分散に対する収束性の比較」(Comparison of con-vergence rate for estimated variance of high dimensional portfolio with long-memory)
2. 13:55-14:35 泉澤 佑 (Izumisawa yu), 白石 博 (Shiraishi Hiroshi)
慶應義塾大学理工学研究科 (Keio University Graduate School of Science and Technology)
慶應義塾大学理工学部数理科学科 (Keio University Department of Mathematics)「Hawkes 過程によるシステミックリスク評価」(Systemic Risk Assessment by Hawkes process)3. 14:35-15:15 藤森洸 (Kou Fujimori), 西山陽一 (Yoichi Nishiyama)
早稲田大学大学院基幹理工学研究科 (Waseda University, Graduate School of Fundamental Scienceand Engineering)
早稲田大学国際教養学部 (School of International Liberal Studies, Waseda University)
The Dantzig selector for diﬀusion processes with high-dimensional parameter
Coﬀee Break 15:15-15:30
Afternoon Session II (in Japanese) 15:30-17:30
Chair Yoichi Nishiyama (School of International Liberal Studies, Waseda University)
4. 15:30-16:10 吉田拓真 (Takuma Yoshida),
鹿児島大学 (Kagoshima University)
Nonparametric smoothing for extremal quantile regression
5. 16:10-16:50 Masanari, MOTOYAMA and Yoshihiko MAESONO,
National Federation of Agricultural Cooperative Associations
Faculty of Mathematics, Kyushu University
On direct kernel estimator of density ratio
6. 16:50-17:30 柿沢佳秀 (Yoshihide Kakizawa),
北海道大学 (Hokkaido University)「対称分布ベース q-MIG カーネル、及び、歪分布ベース q-BS カーネルに基づく、非負データに対す
る密度推定」(Symmetrical-based q-MIG kernel density estimation, and skew-based q-BS kernel

density estimation for nonnegative data)Saturday, 18, November
Morning Session (in Japanese or in English) 9:10-11:50
Chair Takeshi Kato (Faculty of Science and Technology, Sophia University)
7. 9:10-9:50 張元宗 (Chang Yuan-Tsung), 篠崎 信雄 (Shinozaki Nobuo), William E. Strawderman
目白大学 (Mejiro University)
慶應義塾大学 (Keio University)

Rutgers University「未知で等しい分散を持つ p 個の非負な正規母平均の同時推定」(Simultaneous estimation of p

positive normal means with common unknown variance)
8. 9:50-10:30 島谷健一郎 (Ichiro Ken Shimatani),
統計数理研究所 (The Institute of Statistical Mathematics)「方向の piece-wise regression model とその魚群データへの応用」(Circular piece-wise regression

model and its application to movement patterns for school of ﬁshes)
9. 10:30-11:10 石井 晶（Ishii Aki）,
東京理科大学 理工学部 情報科学科 (Department of Information Sciences, Tokyo University of
Science)
Asymptotic properties of classiﬁcation procedures based on eigenstructures in high-dimensional
context
10. 11:10-11:50 Junichi Hirukawa and Sangyeol Lee
Niigata University
Seoul National University
Asymptotic Properties of Mildly Explosive Processes with Locally Stationary DisturbanceLunch 11:50-13:15
Afternoon Session I (in English) 13:15-15:15
Chair Shiraishi Hiroshi (Keio University Department of Mathematics)
11. 13:15-13:55 Daichi Inoue and Takeshi Kato
Graduate School of Science and Technology, Sophia University
Faculty of Science and Technology, Sophia University
Activation Function Selection in Recurrent Neural Network
12. 13:55-14:35 阿部俊弘 (Abe, Toshihiro), 塩濱 敬之 (Takayuki Shiohama), 宮田 庸一 (Yoichi
Miyata)南山大学理工学部 (Nanzan university, Faculty of Science and Engineering)
東京理科大学工学部 (Department of Information and Computer Technology, Tokyo University of

Science)
高崎経済大学経済学部 (Faculty of Economics, Takasaki City University of Economics)
On mode and antimode preserving circular distributions
13. 14:35-15:15 宮田 庸一 (Yoichi Miyata), 塩濱 敬之 (Takayuki Shiohama), 阿部 俊弘 (Toshihiro Abe)
高崎経済大学経済学部 (Faculty of Economics, Takasaki City University of Economics)

東京理科大学工学部 (Department of Information and Computer Technology, Tokyo University of

Science)
南山大学理工学部 (Department of Systems and Mathematical Science, Nanzan University)「超球面上の非対称な分布の有限混合モデルに対する推定について」(On estimating ﬁnite mixturesof skew-rotationally-symmetric distributions)Coﬀee Break 15:15-15:30
Afternoon Session II (Guest speakers session) 15:30-17:30
Chair Junichi Hirukawa (Niigata University)
14. 15:30-16:30 Giuseppe Cavaliere and Iliyan Georgiev
University of Bologna
BOOTSTRAP INFERENCE UNDER RANDOM DISTRIBUTIONAL LIMITS
15. 16:30-17:30 Shiqing Ling, Ruey S. Tsay and Yaxing Yang
Hong Kong University of Science and Technology
University of Chicago
Xiamen University
Testing for Series Correlation and ARCH Eﬀect of High-Dimensional Time Series Data
Conference Dinner 18:00-21:00
ぜんてい 新潟駅南店(Zentei Niigata Ekinanten)
5,500 yen
(TEL: 025-244-8807)

Sunday, 19, November
Morning Session (in Japanese) 9:10-11:50
Chair Abe, Toshihiro (Department of Systems and Mathematical Science, Nanzan University)
16. Shuichi Shinmura,
Seikei University, Prof. Emeritus
Cancer Gene Analysis by Singh et al. Microarray Data
17. 9:50-10:30 佐藤 駿 (Shun Sato), 合田 新平 (Shimpei Goda), 三枝 祐輔 (Yusuke Saigusa),
生亀 清貴 (Kiyotaka Iki), 山本 紘司 (Kouji Yamamoto)
東京理科大学大学院 理工学研究科 (Graduate School of Science and Technology,Tokyo Universityof Science)
東京理科大学大学院 理工学研究科 (Graduate School of Science and Technology,Tokyo Universityof Science)
横浜市立大学 医学部 (Faculty of Medicine,Yokohama City University)) 東京理科大学 理工学部 (Faculty of Science and Technology,Tokyo University of Science)
大阪市立大学大学院 医学研究科 (Graduate School of Medicine,Osaka City University)「二元分割表における独立性および対称性に関するモデルと分解」(Models and decompositions of
independence and symmetry for two-way contingency tables)
18. 10:30-11:10 安藤宗司 (Shuji Ando), 田畑耕治 (Kouji Tahata), 富澤貞男 (Sadao Tomizawa)
東京理科大学工学部 (Faculty of Engineering, Tokyo University of Science)
東京理科大学理工学部 (Faculty of Science and Technology, Tokyo University of Science)
東京理科大学理工学部 (Faculty of Science and Technology, Tokyo University of Science)「正方分割表における対称性の尺度について」(On measures for symmetry in square contingency

tables)19. 11:10-11:50 野間久史（Hisashi Noma）,
統計数理研究所（The Institute of Statistical Mathematics）「ネットワークメタアナリシスによる Comparative Eﬀectiveness Research と高次漸近理論に基づく推測手法」(Comparative eﬀectiveness research using network meta-analysis and eﬀective

inference methods via higher order asymptotics)
Closing 11:50-11:55 Junichi HIRUKAWA (Niigata University)

長期記憶性を持った高次元ポートフォリオの分散に対する収束性の比較

Comparison of convergence rate for estimated varianceof high dimensional portfolio with long-memory
慶應義塾大学理工学研究科
慶應義塾大学理工学部

阿部
白石

文貴
博

1 はじめにポートフォリオ理論とは分散投資を行う際にリスクを回避し, 収益率を高めるには, どう最適化すればよいかを決定する理論であり, 平均-分散基準では, 個々の投資資産の期待収益率 µ と分散共分散行列 Σ の関数で, 最b と標本分散
適分配比率が決定する. 古典的な手法では, µ と Σ の代わりに過去の収益率の標本平均ベクトル µb sam を用いるが, 近年の市場規模の拡大により投資する資産の数が増加し, Σb sam が良い推定量で共分散行列 Σはなくなっている. そこで Fan et al(2008) は収益率にファクター構造を考え, 最小二乗法による Σ の推定量

b LS の方がより収束レートが良い事を示した. しかし, Fan et al(2008) では誤差項が i.i.d. と仮定していたが,Σ
その仮定は現実的ではないと思われる. 本発表では, 誤差項が長期記憶過程に従うと仮定し, その下でのポート
フォリオの収束性について議論する.

2 GMV ポートフォリオとファクターモデル
Yit を資産 i の時刻 t における収益率とする. ここで yt = (Y1t , · · · , Ypt )⊤ は平均ベクトル E(yt ) = µ, 分散
共分散行列 Cov(yt ) = Σ > 0 の分布に従っていると仮定する. さらに配分比率を ξ = (ξ1 , . . . , ξp )⊤ とする. 配分比率はその総和が 1 となるようにする. この時, ポートフォリオの収益率の分散は Var(ξ ⊤ yt ) = ξ ⊤ Σξ となる. この分散が最も小さくなる GMV ポートフォリオの配分比率 ξ ∗ とその時の分散 ξ ∗ ⊤ Σξ ∗ はξ∗ =Σ−1 1p
−1 1
1⊤
p
pΣ

ξ ∗ ⊤ Σξ ∗ =

1
−1 1
1⊤
Σ
p
p

(＊)

これより GMV ポートフォリオの分散を推定する事は Σ−1 を推定することに帰着する.
ファクターモデルは経済工学や金融工学の分野で幅広く用いられている概念である. ファクターモデルの下で収益率 Yit は, 次のように定義される. また次の仮定を置く.

(仮定 1)

p は n と共に増加し p >> K.

(仮定 2)

ft = (f1t , · · · , fKt ) と ϵt = (ϵ1t , · · · , ϵpt )⊤ は無相関.

定義 2.1 (誤差項が長期記憶性を持つファクターモデル)

Yit = b⊤i ft + εiti = 1, · · · , p

ここで f1t , · · · , fKt はファクター (確率変数), bi1 , · · · , biK (i = 1, · · · , p) はファクターローディングと呼ばれる係数である. 特に誤差項 {εit } は長期記憶過程に従い，E(ϵit ) := 0, Var(εit ) := σi2 , Cov(εi0 , εit ) :=σi2 γt ,(t = 1, 2, . . .) とおくと∑∞

t=0 |γt | = ∞ を満たす.

さらに行列形式にまとめて書くと

yt = Bft + ϵt

(2.1)

ここで,B = (b1 , · · · , bp )⊤ である. これから yt の分散共分散行列は

Σ = Cov[yt ] = Cov[Bft ] + Cov[ϵt ]= BCov[ft ]B ⊤ + Σ0

(2.2)

と書けることがわかる. さらに X = (f1 , · · · , fn ), Y = (y1 , · · · , yn ), E = (ϵ1 , · · · , ϵn ) と定める. (2.2) 式か

b LS , Σb BLU , Σb sBLU , Σb F GLS を定める. しかし, 自己共分散行列 Ω を推定し, そこから逆行列をら, Σ の推定量 Σ
考える時, 誤差が大きくなることが知られている. この誤差を小さくするために, Ω−1 を直接推定する.

定義 2.2 (閾値 κ を使った Ω−1 の近似)Ω−1 を修正コレスキー分解によって Ω−1 = T ⊤ D −1 T となるような T , D を定める. ここで T = (tjk )：下三角行列, D = diag(d1 , · · · , dn )：対角行列 である. この時, 閾値 κ を用いて, 次のように T (κ), D(κ) を定める.
0
T (κ) = (tjk (κ)) = tκktjk

if κ + 1 < j ≤ n, 1 ≤ k ≤ j − κ − 1
if κ + 1 ≤ i ≤ n, j − κ ≤ k ≤ j − 1
if otherwise

D(κ) = diag(σ1 , · · · , σκ , · · · , σκ )そしてこれらの値を用いて Ω−1 (κ) = T (κ)⊤ D −1 (κ)T (κ) と定める.
定義 2.3 (FGLSE

実現可能一般化最小二乗推定量)b⊤bb F GLS = BdbΣF GLS Cov(f )BF GLS + Σ0

(2.3)

b F GLS = (X ⊤ Ω̌−1 (κ)X)−1 X ⊤ Ω̌−1 (κ)Y , Σb⊤bbb 0,F GLS = diag(n−1 E
ただし,BF GLS EF GLS ), EF GLS = Y −b F GLS として,Ω̌−1 (κ) は定義 2.2 の Ω に Ω̌ を代入したものである. また κ は Ing et al(2016) に従う.
XB

3 数値結果
紹介した推定量と真の値がどの程度近くなっているかをシミュレーションによって確かめる. シミュレー
ションにはまず 4083 銘柄の対数収益率に対して JASDAQ, 日経平均,TOPIX をファクターとしたファクター
モデルを作る. そしてそこから得られた, ファクターローディングの標本平均 µb と標本分散共分散行列 Σb を計算する. 同時に, ファクターの標本平均 µf と標本分散共分散行列 Σf を計算する. そしてこれらをパラメータとする多変量正規分布からそれぞれファクターローディングとファクターの標本を生成する. そしてこれらの値と ARFIMA(0,d,0) 過程から生成した誤差項の標本を使い, 定義（2.1）で定めたファクターモデルに従って対数収益率の標本を計算する. これらの生成した標本を用いて, 紹介した GMV ポートフォリオの分散の推定量を作成し, パラメータから作成した真の値との差を RMSE によってどの程度乖離しているか評価する. 実際にシミュレーションを行った詳細な結果については当日報告した.

参考文献[1] Ledoit, O., and Wolf, M. (2003). Improved estimation of the covariance matrix of stock returns with anapplication to portfolio selection. Journal of empirical finance, 10(5), 603-621.
[2] Bai, Z., Liu, H., and Wong, W. K. (2009). Enhancement of the applicability of Markowitz’s portfolio
optimization by utilizing random matrix theory. Mathematical Finance, 19(4), 639-667.
[3] Bodnar, T., Parolya, N., and Schmid, W. (2014). Estimation of the global minimum variance portfolio inhigh dimensions. arXiv preprint arXiv:1406.0437.
[4] Goldfarb, D., and Iyengar, G. (2003). Robust portfolio selection problems. Mathematics of operations
research, 28(1), 1-38.
[5] Fan, J., Fan, Y., and Lv, J. (2008). High dimensional covariance matrix estimation using a factor model.
Journal of Econometrics, 147(1), 186-197.
[6] Ing, C. K., Chiou, H. T., and Guo, M. (2016). Estimation of inverse autocovariance matrices for longmemory processes. Bernoulli, 22(3), 1301-1330.
[7] 刈屋武昭, 矢島美寛, 田中勝人, and 竹内啓. (2003). 経済時系列の統計. J 金融 時系列 分析入門, 岩波書店匚, 3,
47-67.

Hawkes 過程によるシステミックリスク評価
Systemic Risk Assessment by Hawkes process
慶應義塾大学理工学研究科 泉澤 佑
慶應義塾大学理工学部
白石 博

1 Introduction
近年、システミックリスクと呼ばれる個別の金融機関のショックが金融システム全体の危機を引き起こすリスクの顕在化が憂慮されている。システミックリスクの重要な要素の一つとして金融危機の伝播・波及というメカニズムがある。金融システム全体の安定性を脅かすようなショックの波及の評価は、システミックリスク研究において重視される。しかし、個々の金融機関でのショックの波及のしやすさ、およびその波及の経路を判断することは困難が生じる。これは波及が取引ネットワークの複雑化、レバレッジなど様々な要因によることに起因し、システミックリスクを評価する際の問題となる。
システミックリスクを定量化するため、金融機関の関連性を取引などでモデル化するなどの研究が行われてきた。[2]では、ショック発生時の悪影響の程度を金融機関同士の取引などでモデル化することでシステミックリスクを定量化している。しかしながらこのようなモデルは悪影響をもたらす直接的、間接的な金融システムの状態・要因全てを捉えられているとは限らない。本研究では Hawkes 過程と呼ばれる点過程モデルを利用することで実経済上に現れる損失発生の波及をモデル化する。このモデルを用いてショックの波及を実経済上に現れる損失の動きから直接、波及経路や波及のしやすさを読み取り、システミックリスクを評価する。

2 Hawkes 過程によるシステミックリスク評価
Hawkes 過程は [3] によって提案された点過程確率モデルである。Hawkes 過程はポアソン型のモデルで、強度パラメー
タに過去のイベントの影響による励起を加えた形で定義 2.1 のように定義される。
定義 2.1. d ∈ N に対して N = (N1 , ..., Nd )T を d 変量 Hawkes 過程とする。このときある d 変量ベクトル η ∈ Rd≥0 、
可測な d 変量正方行列値関数 H = (hi,j )1≤i,j≤d : R≥0 → Rd×d
≥0 に対して d 変量 Hawkes 過程の条件付き強度は[]∫ t
E N ((t, t + δ])|HtN
ΛN (t) := lim=η+H(t − s)N (ds), t ∈ R.
δ↓0
δ
−∞

(2.1)

ここで各 i ∈ {1, ..., d} に対し

(∫ t
−∞

とし、また HtN



)
d ∫ t
∑
H(t − s)N (ds) := hi,j (t − s)Nj (ds) ,
i

j=1

−∞

(

:= σ {ω ∈ Ω : N ((a, b]) = n} , n ∈ Nd0 , a < b ≤ t(hi,j )1≤i,j≤d : R≥0 → Rd×d
≥0 を励起関数行列と呼ぶ。

)

(2.2)

i

とする。このような η ∈ Rd≥0 を基底強度、H =Hawkes 過程の強度は基底強度と過去のイベント発生からの励起の和として表される。H は各成分に関数が入った行列で、例えば Hi2 i1 は i1 でイベントが発生した場合の i2 の強度の励起度合いを表す。
Hawkes 過程の定義は個別のイベント発生の影響が波及し、他でのイベント発生を引き起こすという考えをモデル化したもので、システミックリスク顕在化のメカニズムと非常に合致したものである。金融システムにおけるショックの波及
ネットワークを Hawkes 過程によりモデル化することで、システミックリスクを評価する。Hawkes 過程の波及を評価す
るためその可視化である Hawkes グラフ、そしてそれを元にした二つの指標を導入する。Hawkes グラフとは [1] で提案された Hawkes 過程の波及構造を重み付き有向グラフにより可視化したもので定義 2.2 のように定義される。
定義 2.2. N を基底強度 η 、励起関数行列 H の d 変量 Hawkes 過程とし、K を次のように定義する。

(∫ ∞
K :=
0

)
Hj,i (t)dt
1≤i,j≤dこのとき N に対する Hawkes グラフ G は次の式を満たす頂点集合 VN 、辺集合 EN を持つ重み付き有向グラフとして定義される。

VN := {(j; ηj ) : j ∈ {1, 2, ..., d}; ηj 重み } ,EN := {(i, j; Ki,j ) : (i, j) ∈ VN × VN , Ki,j > 0; Ki,j 重み } .
定義 2.2 において K は Hawkes グラフの隣接行列にあたる。Ki,j >0 であり Hawkes グラフに i から j への辺が存在するとき、つまり j でのイベント発生確率が i の過去の情報に依存するとき i から j の影響が存在することを表し、辺の重みは影響の大きさを表す。頂点の重みは各変量の基底強度であるのでイベント発生の基底となる確率を表す。Hawkes
過程からそれに対する Hawkes グラフを構成することで、波及ネットワークの可視化ができる。
また、Hawkes 過程の特徴を表す指標としてて cascade 係数、feedback 係数を導入する。それぞれ定義 2.3 のように定義される。
定義 2.3. N を subcritical な d 変量 Hawkes 過程とし、K, GN をそれぞれ対応した隣接行列と、それに基づく Hawkes∑
グラフとする。(ei,j ) := g≥0 K g ∈ Rd×d
≥0 として ei,j を定め、これを用いて次のように定義をする。

ηici := ∑d∑d

j=1 ei,j,∑d
l=1 ηlj=1 el,jηj ej,jfj := ∑d, j ∈ {1, 2, ..., d}。
i=1 ηi ei,j{ci }i=1,...,d を cascade 係数、{fj }j=1,...,d を feedback 係数と呼ぶ。
各 i について cascade 係数の大きさは Hawkes 過程全体でのイベント発生と影響の波及全体に対する頂点 i の寄与を表す。対して j での feedback 係数の j が受ける波及のによる影響への j 自身の寄与を表している。cascade 係数、
feedback 係数をそれぞれ金融機関ごとの金融システム全体でのリスクの顕在化への寄与度、金融ネットワークにおける自身へのショックの伝播に対する寄与度と見ることでシステミックリスク評価の指標とすることでができる。
損失波及モデルHawkes グラフはイベントの波及を評価するモデルであった。システミックリスクにおけるショックの伝播による損失の定量的な評価を考える。損失の発生が Hawkes 過程に従うとし、発生した損失に確率分布を考えるマーク付き
Hawkes,Poisson 過程を利用する。定義した損失を、ショックの波及により生じる損失の定量的な評価に用いる。
L(B) := (L1 (B), ..., Ld (B))T , d ∈ N, B ∈ B(R) をある期間 B の経済指標・金融商品の損失を表す確率変数とし、Lを次のように定義する。

  ∑N1∗ (B)  ∑N1 (B)L1 (B)
X1,k
k=1 X1,k
k=1
 

..
..
..
=
L(B) := 

−
.
.
.
.
∗
∑
∑
N
(B)
N
(B)
d
d
Ld (B)
Xd,k
Xd,k


k=1

(2.3)

k=1

∗

ここで N , N をそれぞれ d 次元 Hawkes 過程、d 次元 Poisson 過程に従う確率変数とする。また、Xi,k を N , N ∗ , i, k
について独立な損失を表す確率変数とする。Hawkes 過程はイベント発生の影響による強度パラメータの励起を Poissonに加えたモデルであることを考慮すると、この L は全体の損失から基底強度による、つまり過去の損失の影響を除いた波及によって生じる損失のみを表している。L に対してリスク尺度を適用することによりシステミックリスクにおける波及による損失の定量的な評価とする。

3 数値結果
紹介した指標を金融機関・指標に適用することで、金融システムの脆弱性によるリスクを評価する。経済指標、金融商
品の対数成長率の符号を逆にしたものを負のリターンとする。評価の対象となる期間における負のリターンが適当な閾値を超えた値を損失と定義し、その発生が Hawkes 過程に従い、大きさが一般化パレート分布に従うとする。推定した
Hawkes 過程のパラメータから Hawkes グラフ、cascade 係数、feedback 係数を推定する。または波及による損失 L の期待値をとることでシステミックリスクの定量的な評価を行う。実際の金融市場のデータを解析した詳細な結果については当日報告した。

参考文献[1] Embrechts, Paul, and Matthias Kirchner. “Hawkes graphs.” arXiv preprint arXiv:1601.01879 (2016).
[2] Adrian, Tobias, and Markus K. Brunnermeier. “CoVaR.” The American Economic Review 106.7 (2016): 17051741.
[3] Hawkes, A. G., & Oakes, D. (1974). A cluster process representation of a self-exciting process. Journal of Applied
Probability, 11(3), 493-503.

2

The Dantzig selector for diﬀusion processes withhigh-dimensional parameters
Kou Fujimori

Yoichi Nishiyama

Waseda University
In this talk, statistical models of diﬀusion processes with high-dimensionalparameters were discussed. To estimate high-dimensional parameters, the Dantzig
selector which was proposed by Candés and Tao in 2007 is applied. Then, we
proved two types of consistency of the estimator; one is the lq consistency forevery q ∈ [1, ∞] and the other is the variable selection consistency. Moreover, wetalked about the construction of an asymptotically normal estimator by using
the variable selection consistency of the Dantzig selector.
For example, we considered the following linear model of diﬀusion process:∫ t
Xt = X0 +ΘT ϕ(Xs )ds + σWt ,(1)0where {Wt }t≥0 := {(Wt1 , . . . , Wtp )}t≥0

is a p-dimensional standard Brownian
motion, Θ is a p × p sparse deterministic matrix, and σ = diag(σ1 , . . . , σp ) is ap × p diagonal matrix and ϕ(x) = (ϕ1 (x1 ), . . . , ϕp (xp )) for x = (x1 , . . . , xp ) ∈ Rp
is a smooth Rp -valued function. Note that {Xti }t≥0 for each i = 1, 2, . . . , psatisfies the following equation.
∫ tii
Xt = X0 +ΘTi ϕ(Xs )ds + σi Wti ,0where Θi is an i-th row of matrix Θ. In this talk, we considered the estimation
problem of Θ and σ using the observed data of the process {Xt }t≥0 at n + 1
discrete time points 0 =: tn0 < tn1 < . . . < tnn , where tnk = ktnn /n, k = 0, 1, . . . , n.
Write T0i for the support index set of the true value Θ0i for every i ∈ {1, 2, . . . , p},i.e., T0i = {j : Θ0ij ̸= 0}. Let Si be the number of elements in the index set
T0i . The estimation problem is considered in the following high-dimensional andsparse setting for the true value Θ0 .
p = pn ≫ n,

sup Si =: S ∗ ≪ n,1≤i<∞where S ∗ > 0 is a constant which does not depend on n. It is well known that
we can ignore the influence of Θ when we estimate the diﬀusion coeﬃcients σ.
1

Therefore, we define the estimator σ̂i for σi0 for every i = 1, 2, . . . , pn by the
maximum quasi-likelihood estimator, which satisfies the consistency uniformly
in i. Then, we define the Dantzig selector type estimator for Θ0i by plugging σ̂iin the log-quasi-likelihood as follows:Θ̂n,i := Θ̂i := arg min ∥Θi ∥1 ,Cni := {Θi ∈ Rpn : ∥ψn (Θi )∥∞ ≤ γn },i
Θi ∈Cnwhere γn is a tuning parameter and ψn (Θi ) = ψn (Θi , σi ) is the gradient of
normalized log-quasi-likelihood with respect to Θi . We discussed the asymptotic
property of this estimator. The next theorem states the lq consistency of Θ̂i .
√Theorem 1. Suppose that pn ≫ n, log pn = o( n∆n ), √n∆n = tnn → ∞,22∆n → 0, and n∆n → 0 as n → ∞. Put γn = O(log(1 + pn )/ n∆n ). For global
positive constants K1 , K2 , K3 , the following hold true under some appropriate
conditions.
(i)

(
lim P

n→∞

(ii)

(
lim P

n→∞

sup
1≤i≤pnsup1≤i≤pn(iii)∥Θ̂i − Θ0i ∥1 ≥ K1 S ∗ γn1
∥Θ̂i − Θ0i ∥q ≥ K2 S ∗ q γn

lim P

= 0.

)
q ∈ (1, ∞).

= 0,

(
n→∞

)

)
∥Θ̂i − Θ0i ∥2∞ ≥ K3 γnsup1≤i≤pn= 0.

Then, we talked about the next theorem which states the variable selection
consistency of Θ̂i .
Theorem 2. The estimator T̂ni for T0i defined by1T̂ni = {j : |Θ̂ij | > γn2 }satisfies the following condition under the same assumption as Theorem 1()
lim P T̂ni = T0i for all i ∈ {1, 2, . . . , pn } = 1.
n→∞Finally, we talked about the construction of the asymptotically normal estimator.
(2)

Theorem 3. Define the estimator Θ̂iby the solution to the following equation:
ΘT̂ ic = 0.

ψn (ΘiT̂ i )T̂ i = 0,

(2)

n

n

n

Then, under some appropriate assumptions, it holds for every i that√(2)tnn (Θ̂iT̂ i − Θ0iT i )1{T̂ i =T i } →d N (0, Qi−1),
T i ,T i
n

0

n

0

where QT0i ,T0i is some Si × Si invertible matrix.

2

0

0

Nonparametric smoothing for extremal quantile regression
吉田拓真 (Takuma Yoshida)1
鹿児島大学 (Kagoshima University)
E-mail: yoshida@sci.kagoshima-u.ac.jp
1

1

はじめに

分位点回帰は外れ値に対するロバスト性を有するなど平均回帰にはない利点を持つが, 分位点 τ が
0 や 1 に近い場合 (極値分位点と呼ぶ) はその限りではない. 極値分位点では, 周辺のデータ数が中央に比べ少なくなることから推測が困難となる. その困難さを定式化するためには, 漸近理論と極値理論の融合的な議論を行う必要がある. Chernozhukov (2005), Chernozhukov and Fernández-val(2011), Wang et al. (2012) では線形モデルについて極値分位点回帰の理論結果を与えている. しかし, 中央値が線形構造を有する場合でも, 最大値や最小値に近い分位点が線形構造を有するケー
スは稀である. このような背景から, 講演者は極値分位点回帰のノンパラメトリック推定の議論が重要であると考えた. 本講演では, 極値分位点回帰のノンパラメトリック推定の収束レートと関連する漸近理論結果を報告した.

2

極値分位点回帰

目的変数 Y と説明変数 x に対して, F を Y の分布関数, FY (·|x) を x を与えた下での Y の条件付分布関数とする. また, q(τ ) = F −1 (τ ), qY (τ |x) = FY−1 (τ |x) をそれぞれ Y の 100τ % 分位点, x を与えた下での Y の条件付 100τ % 分位点とする. 極値分位点の理論評価をするためには, 分布の裾の推定の困難さを数学的に扱う必要がある. そこで, データ数 n が大きくなる漸近理論の枠組みで,
τ が 1 に収束するような設定を考える. 今回は n(1 − τ ) → ∞ と n(1 − τ ) → c > 0 の 2 通りの設定
でノンパラメトリック推定量の漸近収束レートを導出する.
まず, 分布の裾について, 以下の仮定をする:

1. y → ∞ の下で, F (y) = y −1/γ L(1 + o(1)) を満たす γ > 0, L > 0 が存在する.
2. ある関数 f , H が存在して, V = Y − f (x) に対して, y → ∞ の下で
1 − FV (y|x) = H(x){1 − F (y)}(1 + o(1))
が成立する. ここで, FV (·|x) は x を与えた下での V の分布関数である.
すると分位点に対して τ → 1 の下で q(τ ) = (1 − τ )−γ L̄(1 + o(1)) とqY (τ |x) = f (x) + H(x)γ q(τ )1 + o(1))
が成立する. この qY (τ |x) を B-スプライン関数で近似することを考える. いま,s(τ |x) = B(x)T b(τ )を p 次 B-スプラインモデルとする. ただし, B(x) は B-スプラインベクトル, b(τ ) はパラメータである. このとき, ある b∗ (τ ) が存在して, s∗ (τ |x) = B(x)T b(τ ) が |s∗ (τ |x) − qY (τ |x)| =1

O(K −m (1 − τ )−γ ) を満たすことが示せる. これは, B-スプラインモデルが持つ性質の簡単な応用である. 具体的には, 同様に f (x) と H(x)γ を近似する B-スプラインモデル sf (x) = H(x)T bf とsH (x) = H(x)T bH を用いて, b∗ (τ ) = bf + q(τ )bH と構成すればよい. つまり, 上記の 2 つの仮定は極値分位点関数をノンパラメトリック推定を行うために極めて本質的な仮定である.

3

ノンパラメトリック推定

以下は報告した研究の主結果である. B-スプラインモデル s(x) = B(x)T b に対して, 推定量をq̃Y (τ |x) = B(x)T b̃(τ ) とする. ただし, b̃(τ ) はn∑ρτ (Yi − s(xi )) + λ∫ b{s(m) (x)}2 dx,ai=1の最小化として求められるものである. ここで, λ > 0 は平滑化パラメータである.
このとき, 簡単な条件と n → ∞, τ → 1, n(1 − τ ) → ∞ の下で[{}2 ]q̃Y (τ |x)E−1= O({(1 − τI )n}−2m/(2m+1) )qY (τI |x)が成立する. 次に, n(1 − τ ) → c > 0 のケースに有用な推定量を構成する. まずは, τ1 > · · · > τk
に対して,
k−1
1 ∑ q̃Y (τj |x)γ̂(x) =k−1
q̃Y (τk |x)
j=1とし, 分位点関数の推定量を(
q̂Y (τ |x) =1 − τI1−τ)γ̂(x)q̃Y (τI |x)とする. ただし, τI は n(1 − τI ) → ∞ を満たすとする.
このとき, 簡単な条件と n → ∞, τ → 1, n(1 − τ ) → c の下で
[{
}2 ]
)
(
{
(
})
2m
2m
q̂Y (τ |x)− 2m+1− 2m+12 1 − τI
E
−1
= O max klog, {(1 − τI )n},qY (τ |x)1−τが成立する.

4

数値例

本講演では, 上記の提案した推定量の精度評価をモンテカルロ・シミュレーションで確認した. また, 実データに対して提案手法を適用した結果を報告した.

References
Chernozhukov,V. (2005). Extremal quantile regression. Ann. Statist. 33 806–839.
Chernozhukov,V. and Fernández-val. (2011). Inference for extremal conditional quantile models,with an application to market and birthweight risks. The Review of Economic Studies. 78 559–
589.
Wang,H.J., Li,D. and He,X. (2012). Estimation of high dimensional conditional quantiles forheavy-tailed distributions. J. Amer. Statist. Assoc. 107 1453–1464.
2

On direct kernel estimator of density ratio
全国農業協同組合連合会
九州大学大学院数理学研究院

本山真誠
前園宜彦

本講演では密度関数の比のカーネル型推定量について議論する．X1 , . . . , Xm と Y1 , . . .
, Yn をそれぞれ分布関数 F (x)，G(x) 及び密度関数 f (x)，g(x) を持つ母集団からの無作為標本とする．このとき，点 x0 ∈ R での密度関数比 f (x0 )/g(x0 ) の推定を考察する．ただし g(x0 ) ̸= 0とする．簡単のために m = n と仮定する．自然な推定量としては，分子・分母にそれぞれのカー
ネル型密度関数推定量を代入するものが提案されている．カーネル関数 K を∫K(u)du = 1
を満たし，K(−u) = K(u) とする．このときそれぞれの密度関数のカーネル型推定量は

fb(x0 ) =n
1 ∑K
nh i=1()

x0 − Xi,hgb(x0 ) =n
1 ∑K
nh i=1(x 0 − Yih)で与えられる．ただし h はバンド幅で h → 0, nh → ∞ と仮定する．ここでは簡単のために同じ
カーネル関数を使用する．これを利用すると密度比の推定量として
bf (x0 )
fd(x0 )  b,g (x0 )=
g(x0 )0,

gb(x0 ) ̸= 0
gb(x0 ) = 0

が利用できる．
これに対して密度関数の比 f (x)/g(x) は点 G(x) における G(X1 ) の密度関数となることを利用して Ćwik & Mielniczuk (1989) は直接型推定量
n
1 ∑
fg
(x0 )
K
=
g(x0 )
nh i=1(Gn (x0 ) − Gn (Xi )h)を提案している．ここで Gn (·) は経験分布関数

Gn (y) =n
1∑
I(Yi ≤ y)n i=1

で，I(·) は定義関数である．この推定量の漸近平均二乗誤差 (AM SE) は h = O(n−d )(1/6≤ d ≤ 1/4) のとき
{

AM SE

fg
(x0 )
g(x0 )}={ ′′f (x0 )

2g(x0 )
1
+
nh

{

}

−

2
f g ′′(x)µ22 h40
2g 2

f (x0 )
f 2 (x0 )+g 2 (x0 ) g 3 (x0 )

}∫

1

(

1
K (u)du + O(h ) + o
nh
−1
2

)

5

で与えられる．Gn は階段関数であるから，上記の直接型推定量は離散的な値を取るものになる．
カーネル型推定量は滑らかな推測結果が得られることが重要な利点である．この弱点を補うために，経験分布関数をカーネル型分布関数推定量で置き換えた
⋆

n
fg
(x0 )
1 ∑
=
K
g(x0 )
nh i=1(b 0 ) − G(X
b i)
G(x
h

)

を議論する．ここで

1
b
G(x)
=

n
∑

n i=1

(

W

)

∫ t

x − Yi,h

W (t) =−∞

K(u)duである．この推定量は滑らかなものであり，適当な条件の下で下記の漸近平均二乗誤差を得ることができる．
定理 1 (i) g ′ (u) は存在し |g ′ (u)| ≤ M (M > 0), u ∈ R. (ii) K は対称カーネルでサポートは∫∫[−1, 1] で， u2 K(u)du < ∞, K ′′ (u)du = 0, K (4) (u)| ≤ M とする．もし h = O(n−d ) (1/6 <d < 1/4) ならば提案する推定量の AM SE は
{

AM SE

fg
(x0 )
g(x0 )

⋆}

{ ′′

=

}

2
r (G(x0 )) f g ′′− 2 (x0 ) µ22 h4
2
2g

1
+
nh

{

f (x0 ) f 2 (x0 )+ 3
g(x0 )
g (x0 )となる．ただしr(u) =}∫

1(

1
K (u)du + O(h ) + o
nh
−1
2

)

5

f (G−1 (u))
g(G−1 (u))

である．
また次の漸近正規性も成り立つ．
定理 2 定理 1 の仮定が成り立ち，もし h = n−1/5 ならば

√

[

⋆

]

fg
(x0 )
f (x0 ) Lnh−−→ N (µ⋆ , σ⋆2 )g(x0 )g(x0 )

が成り立つ．ただし
}

{ ′′

µ⋆ =
{

σ⋆2

=

r (G(x0 )) f g ′′− 2 (x0 ) µ2 ,
22g
f (x0 ) f 2 (x0 )+ 3
g(x0 )
g (x0 )}∫

1−1

K 2 (u)du

である．

参考文献
[1] Chen, S. M., Hsu, Y. S., and Liaw, J. T. (2009). On kernel estimators of density ratio.
Statistics, 43, 463-479.
[2] Ćwik, J., and Mielniczuk, J. (1989). Estimating density ratio with application to discriminant analysis. Communications in Statistics-Theory and Methods, 18, 3057-3069.
[3] Nadaraya, E.A. (1964). On estimating regression, Theory of Probability and Its Applications,10, 186-196.
[4] Watson, G.S. (1964). Smooth regression analysis, Sankhyā Ser. A, 26, 359-372.

対称分布ベース q-MIG カーネル, 及び, 歪分布ベース q-BS カーネルに基づく, 非負データに対する密度推定
柿沢 佳秀 (北大経済)1. はじめに X1 , . . . , Xn ∼ iid f とする. Rosenblatt (1956) の (標準的な) カーネルi
) の漸近的性質はよく知られている (Silverman (1986),
密度推定量 (nh)−1 ni=1 K( x−XhWand and Jones (1995)) が, 推定密度 f の台が R でないならば一般に境界バイアスがあり (一致性をもたない), この対処法も様々に議論された (リノーマライゼーション法,
リフレクション法, 一般化ジャックナイフ法は Jones (1993) により詳しく解説された;
なお, このカテゴリーには “局所多項式法による密度推定”も含まれる).
一方, 推定密度 f が区間 [0, 1], 半無限区間 [0, ∞) = R+ の台を持つ場合, Chen (1999,
2000) がベータ, ガンマカーネル密度推定量を提案して以来, 『非対称カーネル法』が脚光を浴びている. 平滑化パラメータを b = bn > 0 とし, この方針から現在までに議論されている, (点 x ∈ R+ での) 密度推定量 n−1 ni=1 kx,b (Xi ) を３つに分類しておく1.
(1) ガンマ, 逆ガンマカーネル密度推定量を含む Amoroso カーネル密度推定量族(Igarashi and Kakizawa (2017)).
(2) バーンバウム・サンダース (BS), 逆ガウス (IG), 相反逆ガウス (RIG) カーネル
密度推定量を含む, IG と RIG との -混合 (MIG) カーネル密度推定量族 (Igarashi andKakizawa (2014)), 及び, 対数正規 (LN), 加重 LN カーネル密度推定量 (Igarashi (2016));
なお, MIG と LN は “q-MIG”として統一的に扱える (すなわち, LN は q → 0 に相当).
(3) インデックスが ν = ±1/2 の一般化逆ガウス (GIG) カーネル密度推定量 (なお,ν = −1/2(IG), 1/2(RIG) を (2) に分類した; Igarashi and Kakizawa (2014) は GIG の代わりの別名 “変形ベッセル”を採用).
ところで, (2) のカーネルが標準正規密度 φ を元に構成されており, 本報告では次の2 点から, 非負データに対する密度推定量族を拡張した2:(A) 標準正規密度 φ(u) を対称密度 Cg g(u2) に変更することから, 対称分布をベースにした q-MIG カーネル密度推定量族 (Kakizawa (2018)) への拡張 (なお, q → 0 は対数対称カーネル密度推定量族になる).
(B) (A) の部分族である “BS 型”を, ジョンソンシステム (Johnson (1949)) の分布の変換法のアイディアにより拡張 (なお, この拡張は IG 型, RIG 型では不可能である).
1

このような形式は “重み付き推定量” (例えば, Siliverman (1986, Subsection 2.9: General weightfunction estimators)) であるが, Chen (1999,2000) 以降は “非対称カーネル推定量”がキーワードである.
注意. いくつかの推定量は, オリジナルの提案 (Jin and Kawczak (2003), Scaillet (2004), Koul and Song(2013), Marchant et al. (2013), Saulo et al. (2013) など) において(i) fb (0) = 0 の制約 (明らかに, f (0) > 0 の場合には適さない)(ii) O(1) の境界バイアスの問題(iii) 積分分散の主要項の存在性に関わる問題 (その存在を仮定すると, f (0) > 0 なる R+ 上の有界な密度関数が排除される)などがあり, いずれかの問題があるものについては, それらが回避された “再定義”も議論されている.
2
報告では非負データに対するノンパラメトリック密度推定問題への応用を念頭にしているが, 関連する分布論の検討でもあることに注意したい.

2. 分布論 (A) Dı́az-Garcı́a and Leiva-Sánchez (2005) は BS 密度において, φ(u) を対称
密度 Cg g(u2) に替え, GBS 分布 BSg (α, β) を定義した (報告者は, 対称分布をベースにした BS 分布と名付けた). 一方, Yang (2006) は, dual transformationaq (t) =⎧ qt − t−q
⎪
⎪
⎨
,

2q
⎪
⎪
⎩ log t ,q = 0 ,(strictly increasing from (0, ∞) onto R)

q=0(see Slymen and Lachenbruch (1984)) に基づいて, 古典的な (正規分布をベースにした)BS 密度を拡張することから, q → 0 として LN 密度へ繋いでいる. 両者のアイディアを組み合わせることで, “対称分布をベースにした q-BS 密度”が得られる:
kgBSq (s; α, β) = Cg g {aq (s/β)}2α2

Aq (s/β),αβ

s ∈ R+ .

[BSgq ]ここに, 1/Cg = R g(u2) du = R+ u−1/2 g(u) du, Aq (t) = (d/dt)aq (t).
この [BSgq ] を経由して, 以下のような [IGgq ], [RIGgq ] を導入するには, 積分の等式
R+

g((sq − s−q )2 /(2qα)2 )sq−1 ds =R+g((tq − t−q )2 /(2qα)2 )t−(q+1) dt (s = 1/t)が役立つであろう. すなわち,

kgIGq (s; α, β) = Cg g
kgRIGq (s; α, β) = Cg g {aq (s/β)}2α2
 {aq (s/β)}2α2

1  β q+1,
αβ s[IGgq ]1  s q−1  s 2q IGq=
kg (s; α, β) .
αβ ββ[RIGgq ]なお, g(y) = e−y/2 = gN (y) (say), q = 1/2 の場合が, Tweedie (1957; (1a)–(1d)) による
古典的な (正規分布をベースにした) IG, RIG 密度であり, 古典的な (正規分布をベースにした) BS 密度は, これらの 1/2-混合としても得られる (Jørgensen et al. (1991)).
(B) R 上の密度 h をもつような確率変数 Z(h) ∼ h に対し, ξ(h) = βa−1 (α(Z(h) + θ)) の従う分布として『BS 型分布』を定義すると, その密度はBSk(h)(s; α, β, θ) = ha(s/β)A(s/β)−θ,ααβs ∈ R+になる (この q バージョンも同様に定義). 本報告では, h として, Azzalini (1985) 型のskew 密度, 及び, ツーピース型の epsilon-skew 密度の２種類を適用した (現在進行形).

3. 非負データに対する密度推定 対称分布ベースの q-MIG カーネルkgMIG q (s; α, β) = (1 − )kgIGq (s; α, β) + kgRIGq (s; α, β) ,BSq ∈ [0, 1]

及び, 歪分布ベースの q-BS 型のカーネル k(h) (s; α, β, θ) に付随するパラメータ (α, β)を推定位置 x ∈ R+ に依存させて, ノンパラメトリックな (可変的な) 対称分布ベース
q-MIG カーネル, 及び, 歪分布ベース q-BS カーネルに基づいた密度推定量族を提案し,それらの漸近的性質と数値実験の結果を報告した.

未知で等しい分散を持つｐ個の非負な正規母平均の同時推定 (Simultaneous Estimationof p Positive Normal Means with Common Unknown Variance)
目白大学

張

元宗,

慶應義塾大学 篠崎 信雄 ,Rutgers 大学

William E. Strawderman

1．はじめに X1 , . . . , Xp を互いに独立に正規分布 N (θi , σ 2 ), i = 1, . . . , p にしたがう確率変数とし、θi ≥ 0 とする。等分散 σ 2 が未知で、標準化２乗誤差損失関数を基準としたとき、母平均 θ = (θ1 , · · · , θp )′ の同時推定問
題を考える。
分散が既知で、1 次元の場合、Katz(1961) は θ の事前分布を非負な区間の上の一様分布とするとき、θ の一
般化 Bayes 推定量
δ(X) = X + σψ(X/σ)(1.1)
を提案し、δ(X) は θ のミニマックス推定量であり、許容的な推定量でもあることを示した。ここで、

e−x /2,e−t2 /2 dt−∞
2

ψ(x) = ∫ x

ψ(0) =√

2/π,である。全ての −∞ < x < ∞ に対して、δ(x), ψ(x) の性質は

a) x + ψ(x) ≥ 0,b) ψ ′ (x) = −ψ(x)(x + ψ(x)) ≤ 0,c) x + ψ(x) は x の増加関数である。.
である。
p 次元の場合に対しても δ(X) = (δ(X1 ), . . . , δ(Xp )) は θ の一般化 Bayes 推定量である。しかし、θ を同時推定するとき、標準化２乗誤差損失関数の下で、張 (1982) は δ(X) を改良する縮小推定量
(
)
cσ 2
c
δi (X) = 1 −δ(Xi )(1.2)||δ(X)||2を提案し、 δ c (X) = (δ1c∑(X), . . . , δpc (X)) は δ(X) よりよくなるための十分条件 p ≥ 3, 0 ≤ c ≤ 2(p − 2) を与えp2
た。ここで ||δ(X)|| = i=1 δ 2 (Xi ) である。
本研究の目的は分散 σ 2 が未知の場合、分散が既知の場合の推定量 (1.1) に分散の推定量を代入して得られる推定量√√
δ(Xi , S) = Xi + a Sψ(Xi /(a S))に基づく縮小推定量(
δic (X, S) =

1−

)
cS
δ(Xi , S),||δ(X, S)||2(1.3)

を提案し、δ(Xi , S), i = 1, . . . , p を改良するための十分条件を与える。ここで、S ∼ σ 2 χ2n は Xi と互いに独立とする。さらに、δic (X, S) を含む２つの Baranchik タイプ縮小推定量のクラスを提案する。証明するため、下記の予備定理が有用である。
予備定理 2.1 全ての −∞ < X < ∞ に対して、

a) Sψ(X/S) は S の増加関数である。
b) δ(X, S) は S の増加関数である。(よって、 ||δ(X, S)||2 は S の増加関数である).
√√
c) −Sψ ′ (X/(a S)) は S の増加関数である。
√√√√
d) S −1 (X + Sψ(X/ S)) は S の減少関数である。(よって、 ||δ(X, S)||2 /S 2 は S の減少関数である。)
が成立する。
予備定理 2.2 (χ2 identity)

S ∼ σ 2 χ2n とする。
E[h(S)Sσ 2 ] =E[S 2 h(S)] − 2σ 2 E[S 2 h′ (S)].
n+2

が成立する.。
予備定理 2.3 g(s) は s0 で負の値から正の値に変わって、 h(s) は非減少関数であるとするとE[g(S)h(S)] ≥ E[g(S)]h(s0 ).
が成立する。
２．結果
定理 2.1. 標準化 2 乗誤差損失のもとで、δ c (X, S) が δ(X, S) を改良するための十分条件は

p−2 2
1
, a ≥
, n ≥ 5.
n+2n−4
上の結果を下記のような２つの Baranchik タイプ推定量のクラスに拡張する。
タイプ１：
(
))2
p (
∑
Xi
X
√ + ψ √iF = ||δ(X, S)||2 /S = a2
とし、
a S
a S
i=1
(
)
r(F/S)
δiB (X, S) = 1 −δ(Xi , S)
Fとする。次の定理が得られる。
p ≥ 3, 0 < c < 2定理 2.2. δ B (X, S) が δ(X, S) を改良するための十分条件は下記のようである。

i) r(F/S) は F/S の単調非減少関数ii) r(F/S)/(F/S) は F/S の単調非増加関数
p−2iii) 0 ≤ r(F/S) ≤ 2 n+2, a2 ≥ 1/(n − 4), n ≥ 5.

標本の大きさに対する条件を改良するため、次のような推定量のクラスを提案する。
タイプ 2：
{
(
)
}2
p
∑
√
Xi
Xi
2
F̃ = amax √ + ψ √
, 2/π
≥F
a S
a S
i=1
を定義し、

(
δiB F̃ (X, S) =1−)r(F̃ )
δ(Xi , S)F̃

とする。次の定理が得られる。
定理 2.3. δ B F̃ (X, S) が δ(X, S) を改良するための十分条件は下記のようである。

i′ ) r(F̃ ) は F̃ の単調非減少関数
ii′ ) r(F̃ )/F̃ は F̃ の単調非増加関数
p−2iii′ ) 0 ≤ r(F̃ ) ≤ 2 n+2, a2 ≥ 1/(n − 2), n ≥ 3.

参考文献
(1) Baranchik, A. J., 1970. A Family of Minimax Estimation of the Mean of a Multivariate Normal Problems.
Annals of Math. Stat., Vol.41, pp. 642-645.
(2) Chang Y.-T., 1982. Stein-Type Estimators for Parameters in Truncated Space. Keio Science and Technology Reports, Vol. 35, No.11, pp. 185-193.
(3) Birnbaum, Z.W., 1942. An Inequality for Mill’s Ratio. Ann. Math. Stat. Vol. 13, pp. 245-246.
(4) Gordon, R.D., 1941. Values of Mill’s Ratio of Area to Bounding Ordinate and of the Normal Probability
Integral for Large Values of the Argument, Annals of Math. Stat., Vol.12, pp. 364-366.
(5) Katz, M. W., 1961. Admissible and Minimax of Parameter in Truncated Spaces. Annals of Math. Stat.,Vol. 32, pp 136-1142.
(6) Sampford, M. R., 1953. Some Inequalities on Mill’s Ratio and Related Functions . Annals of Math. Stat.,Vol. 24, pp 130-132.
(7) Stein, C. M., 1981. Estimation of the Mean of Multivariate Normal Distribution. The Annal of Statistics,
Vol. 9, No.6, pp. 1135-1151.

方向の piece-wise regression model とその魚群データへの応用
島谷健一郎（統計数理研究所）

Circular piece-wise regression model and its application to movement patterns forschool of fishes
Ichiro Ken Shimatani (The Institute of Statistical Mathematics)
When we investigate relationships between two variables, drawing a scattering diagram is acommon first step, followed by an application of a regression model. When variables are both
circular, this first step does not always give us good visualization, which facilitates theimportance of the second step. Currently, couple of regression models are proposed and havebeen applied to circular-circular data. Among them, Downs and Mardia (2002) proposed thefollowing model;
y = α + β + 2 tan −1{w tan(( x − α ) / 2)} + e ,(1)

where x and y are circular variables (–π ≤ x, y < π), w (–1 ≤ w ≤ 1) is the regression coefficient,
β is the intercept (practically, rotation), and α is the “focal direction” (–π ≤ α, β < π) (because
directions do not have a specific origin, α is included for playing a role of an origin). e is arandom variable from a circular probability distribution such as the von Mises and the wrapped
Cauchy distribution
This model is an example of a generalized linear model (GLM), where tan((x − α)/2) is used
as a link function for transforming circular variables to linear variables. In fact, substituting X =tan((x − α)/2) and Y = tan((y − α)/2) into the linear regression equation Y = wX, we obtainedtan((y − α)/2) = wtan((x − α)/2), which is equivalent to equation (1). A difference from common
GLMs is that the circular–circular regression uses the link function twice because both
independent and dependent variables are circular. Consequently, dependent variables tend to be
attracted to α + β, and the attraction becomes stronger if w is close to 0. If w = 1, (1) reduces toy = x + β + e, (2)which expresses a rotation of angle β.
Occasionally, circular-circular data exhibit more complex, irregular patterns. Here, we focus onwhen a basic pattern is just a rotation, but (x, y)s in some parts do not follow it (even including
stochastic uncertainty), and show concentration on each specific angles.
One method for producing a well-fitting model is to connect several regression models
continuously at the border (piece-wise regression). For circular-circular regression, however,
because of its periodicity, it is not so straightforward.
In this study, I propose the circular-circular regrssion model that is defined as a weighted sumof model (1). The weights are given by functions on a circle that have high peaks in eachsepcific parts.
The expected direction can be expressed by their weighted average over circular variables(the argument of a weighted sum over unit vectors)Jy = arg(∑ g j ( x) f j ( x)) + e ,

(3)

j =0

where f j ( x) = α j + β j + 2a tan(w j tan( X −2α )) (j = 1, …, J), g j (x) is a probability density function ofj

the von Mises or wrapped Cauchy distribution (subtracted by the minimum so that min(g j (x)) =0); high around α j and low (almost = 0) when distant from α j . In addition, we prepare a basicrotation; f 0 (x) = x − β 0 (rotation) and g 0 (x) = c 0 (constant). Dependent variable y is deformed
toward α j + β j when x is around α j , and the degree of deformation is strong if w j is close to 0. If
x is distant from any of α j , it is closed to the rotation of β 0 .
The unknown parameters are (C x , C y , ρ in Eq. (2), c 0 , {α j , β j , ρ j , w j }).
The model was formulated under a Bayesian model (non-informative prior distributions wereused), and we estimated the unknown parameters by Markov chain Monte Carlo (MCMC)
method.
The integer J, the number of components, expresses how many abnormal parts there were in a
swarm. Because model (3) is a kind of “mixture model”, and it is well known that for such a“singular model”, commonly used information criterion such as AIC (Akaike information
criterion) is not applicable. This study used the recently proposed information criterion, widely
applicable Bayesian information criterion (WBIC, Watanabe (2013)), which is for Bayesianmodels and applicable to singular models.
What is important is that a model and parameter estimates should have interpretations
associating with given data and the background. The proposed model is for analyzing
spatio-temporal data about a rotating school of fishes in a aquarium, which was produced anddigitalized by Terayama et al. (2015). I have demonstrated how the model and its parameter
estimates are interpreted from zoological viewpoints.
T.D. Downs and K.V. Mardia. (2002). Circular Regression. Biometrika, Vol. 89, pp. 683-697.
K. Terayama, H. Hioki, and M. Sakagami. (2015). Int. J. Semantic Computing, Vol 9,
pp.143-168.
S. Watanabe. (2013). A Widely Applicable Bayesian Information Criterion. J. Machine Learning
Research, Vol 14, pp. 867-897.

Asymptotic properties of classification procedures based on
eigenstructures in high-dimensional context
Aki Ishii
Department of Information Sciences, Tokyo University of Science, Chiba, Japan

1 Introduction
Suppose we have two classes πi , i = 1, 2, and define independent p × ni data matrices, X i =[xi1 , ..., xini ], i = 1, 2, from πi , i = 1, 2, where xij , j = 1, ..., ni , are independent and identically distributed (i.i.d.) as a p-dimensional distribution with a mean vector µi and covariance matrix Σi (≥ O).
We assume ni ≥ 3, i = 1, 2. The eigen-decomposition of Σi is given by
Σi = H i Λi H Ti =p∑λs(i) hs(i) hTs(i) ,s=1where Λi = diag(λ1(i) , ..., λp(i) ) having λ1(i) ≥ · · · ≥ λp(i) (≥ 0) and H i = [h1(i) , ..., hp(i) ] is an
1/2

orthogonal matrix of the corresponding eigenvectors. Let X i − [µi , ..., µi ] = H i Λi Z i for i = 1, 2.
Then, Z i is a p × ni sphered data matrix from a distribution with the zero mean and identity covariancematrix. Let Z i = [z 1(i) , ..., z p(i) ]T and z j(i) = (zj1(i) , ..., zjni (i) )T , j = 1, ..., p, for i = 1, 2. Notethat E(zjk(i) zj ′ k(i) ) = 0 (j ̸= j ′ ) and Var(z j(i) ) = I ni , where I ni denotes the ni -dimensional identity
matrix. Also, note that if X i is Gaussian, zjk(i) s are i.i.d. as the standard normal distribution, N (0, 1).
We assume that the fourth moments of each variable in Z i are uniformly bounded
for i = 1, 2. Let−1 ∑niTz oj(i) = z j(i) − (z̄j(i) , ..., z̄j(i) ) , j = 1, ..., p; i = 1, 2, where z̄j(i) = ni
k=1 zjk(i) . We assume
that P (limp→∞ ||z o1(i) || ̸= 0) = 1 for i = 1, 2, where || · || denotes the Euclidean norm.
Let x0 be an observation vector of an individual∑
belonging to πi (i = 1, 2).∑We assume x0 and xij sii(xij − xini )(xij −are independent. We estimate µi and Σi by xini = nj=1
xij /ni and S ini = nj=1
Txini ) /(ni − 1). A typical classification rule is that one classifies an individual into π1 if{ det(S ) }2n2(x0 − x1n1 )T S −1
(x
−
x
)
−
log
0
1n
1
1n1
det(S 1n1 )< (x0 − x2n2 )T S −1
2n2 (x0 − x2n2 ),(1.1)

and into π2 otherwise. However, the inverse matrix of S ini does not exist in the HDLSS context (p > ni ).
Aoshima and Yata (2011, 2015) considered substituting {tr(S ini )/p}I p for S ini by using the difference of a geometric representation of HDLSS data from each πi . Aoshima and Yata (2015) gave themisclassification rate adjusted classifier for multiclass, high-dimensional data whose misclassification
rates are no more than specified thresholds under the following condition for eigenvalues:
λ21(i)→ 0 as p → ∞ for i = 1, 2.
(1.2)
tr(Σ2i )Recently, Aoshima and Yata (2017) considered the “strongly spiked eigenvalue (SSE) model” as follows:{ λ2 }
1(i)lim inf> 0 for i = 1 or 2.
(1.3)
p→∞
tr(Σ2i )On the other hand, Aoshima and Yata (2017) called (1.2) the “non-strongly spiked eigenvalue (NSSE)
model”.

2

A new classification under the SSE model

As for the two-sample test, Ishii (2017a, 2017b) gave new test procedures for high dimensional data
under the SSE model. In this talk, we considered the distance-based classifier under one of the SSE
models. Aoshima and Yata (2014) considered a classification rule given by using the identity matrix I p
instead of S ini in (1.1). They showed the asymptotic normality of the classifier and provide a sample
size determination so as to control misclassification rates being no more than a prespecified value. They
further developed the classifier to multiclass classification. In order to remove the bias, we considered adata transformation and gave the following new classification rule. One classifies an individual into π1 if(∑tr(S ini A)
x1n1 + x2n2 )TA(x2n2 − x1n1 ) +
(−1)i
<0
2
2ni
2

x0 −

(2.1)

i=1

and into π2 otherwise. Here, the second term of (2.1) is a bias-correction term. If h1(1) = h1(2) (=h1 , say), we gave A = I p − h1 hT1 . We used the noise-reduction (NR) methodology given by Yata andAoshima (2012) to estimate the first eigenvalues and eigenvectors. By using the asymptotic propertiesof the NR estimator given by Ishii, Yata and Aoshima (2016), we gave the asymptotic properties of our
new classification rule. We also compared classification procedures by using computer simulations.

References[1] Aoshima, M. and Yata, K. (2011). Two-Stage Procedures for High-Dimensional Data. Sequantial
Analysis (Editor’s special invited paper) 30: 356-399.
[2] Aoshima, M., Yata, K. (2014). A distance-based, misclassification rate adjusted classifier for multiclass, high-dimensional data. Annals of the Institute of Statistical Mathematics, 66, 983–1010.
[3] Aoshima, M. and Yata, K. (2015). Asymptotic normality for inference on multisample, highdimensional mean vectors under mild conditions. Methodology and Computing in Applied Probability 17: 419-439.
[4] Aoshima, M. and Yata, K. (2017). Two-sample tests for high-dimension, strongly spiked eigenvalue
models. Statistica Sinica, to appear.
[5] Ishii, A., Yata, K., and Aoshima, M. (2016). Asymptotic properties of the first principal component
and equality tests of covariance matrices in high-dimension, low-sample-size context. Journal of
Statistical Planning and Inference, 170: 186-199.
[6] Ishii, A. (2017a). High-dimensional two-sample tests for non-gaussian strongly spiked eigenvalue
models, to appear in Journal of the Japan Statistical Society.
[7] Ishii, A. (2017b). A two-sample test for high-dimension, low-sample-size data under the strongly
spiked eigenvalue model, Hiroshima Mathematical Journal, 47, in print.
[8] 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.

                                         







 

 !    " #$%&'    #$%$'( ) 

#'

¾

    +          ,    

½+      

 





½





½ '#

½

¾

*  ½½' -#   '    





#

 

      !( ./         ,                          #'-    #$%$'

            ,                ,(  !( 0

Æ     

1 

 

   ,     

 #2334 '          # ½ '( 5      *     

#'( "   3      !                                (             0

1  #2334 '                    +   

#2'

2



.      0 1    #  '  

        ¼

 #2334 '      ¼

# '( . 

   (  




6  #2334'

                 ,         



 

( 1 





0 #233&'     

       





1 

 #232'         7   ( 8 ( #234'            ,       ,                  (              ( .                             (                     0    #0 #$9%''           7     
$$4'       

   

 :   #$$9  $$9  $$9     (                    ( :   #232'      
( 5          ( .       







     ;       ;        

 +   " #$&3'    " #$$$'(

            

 #232'( <  





                 ;                 #233$'  0  ( #23'( <         

 0 =         : 7< 



7     ! 

     (         5            













     



 ¼ 










# '







¼

 





 

½































 - 























 

 #' 

 # '  
¼½
 ¼












   # '
¼½
 ¼







¼





¼

½




½

      7    >3?
  >3 ?
   1 Æ         
¼

   3 ¾ 

  3



 



        1 #'  - 

 

½





½









    7

-







,      (Æ-"      *     
#



 -

#@'



   











     





½

¾


½
'

½


#



'

½


 


 

 ½



¾

  





#A'





-




 



½





        



2

 

½

¾

-







½

 





 ´½ ½µ #3' #'  ´¼ ½µ  







     

















     


 






½







   ( "        (        ,    # ;'   B#' 

. 

 ½  ½ ¾
 ½  ½½

 ½ 

 


 

 

 

2









Æ-Activation function selection in Recurrent Neural Network
Daichi Inoue∗

1

Takeshi Kato†Introduction

The problem of solving the inverse of a time-varying matrix arises in many fields, science, engineering andbusiness, for instance. Even when the matrix is time-invariant, however, it is well known that the actual
computation is not easy if the matrix size is large or sparse. To say nothing of time-varying matrices.
A new approach by applying recurrent neural network was proposed by Zhang et al [1] to computethe inverse of a time-varying matrix, and a series of research works has been made in several papers suchas [3], [4], and [5]. Recurrent neural network (RNN) is one type of neural networks. The advantage of
RNN is that it can handle continuous information such as a sentence, so that it has originally applied tonatural language processing and shown high performance. The research works by Zhang and others, [1],[2], [3], [4], and [5], make use of what it calls ”super-sigmoid function” as a good activation function in
RNN, but the mathematically rigorous argument has been skipped in those papers. Then the aim of ourresearch work is to clarify how to select good activation functions of RNN in theory and propose a newactivation function with higher performance than the super-sigmoid function. The argument is made inview of the convergence of error terms to zero and robustness caused by the diﬀerentiation error and themodel-implementation. Numerical examples are also given in addition to theories.

2

The summary of previous studies

2.1

Problem formulation

We consider computing the inverse of a smooth time-varying matrix A(t) ∈ Rn×n . In other words, wewill find X(t) ∈ Rn×n satisfying the equation
A(t)X(t) − I = 0, t ∈ [0, ∞)(1)

where I ∈ Rn×n is the identity matrix.

2.2

General RNN Model

We define a matrix-valued error function by
E(X(t), t) = (eij (t)) := A(t)X(t) − I
to monitor the matrix-inversion process. The error function derivative Ė(X(t), t) with respect to t
should be made such that every entry eij (t) of E(X(t), t) converges to zero. Specifically, E(X(t), t) canbe described in the general formdE(X(t), t)= −ΓF(E(X(t), t))(2)dt
where F(·) : Rn×n → Rn×n denotes a matrix mapping of neural networks, and the matrix Γ consistingof positive elements scales the convergence rate of the solution. Each element of F = (fij ) is called an
activation function. The RNN design formula (2) leads to the following implicit dynamic equation of thegeneralized neural model
˙
A(t)Ẋ(t) = −A(t)X(t)− ΓF(A(t)X(t) − I)(3)

where X(t), starting from an initial condition X(0) = X0 ∈ Rn×n , is the activation state matrix corresponding to the theoretical solution X ∗ (t) of (1).
∗ Graduate School of Science and Technology, Sophia University (Master’s Program)† Faculty of Science and Technology, Sophia University

3

Comparison of activation functions

We assume below that Γ = γI for some γ > 0 and fij ≡ f in F = (fij ) for simplicity. Zhang et al. [2]discusses the performance of several activation functions in RNN, and shows that the following powersigmoid activation function f1 is the best among them in the sense of convergence of errors to zero andthe robustness: p: |u| ≥ 1
 u
f1 (u) =eξ + e−ξ ξtanh(ξu) : |u| < 1e − e−ξwhere ξ > 1/2 is a real number and p ≥ 3 is an odd integer.
We propose a new activation function1sgn(u)(ep|u| − 1)
 pe −1
f2 (u) =eξ + e−ξ ξtanh(ξu)e − e−ξ:|u| ≥ 1

: |u| < 1with ξ > 1 and an odd integer p ≥ 3.
Solving the diﬀerential equation (2) for f = f1 , f2 and noting the diﬀerence between the power function
and the exponential function in |u| ≥ 1, we can show that that f2 is slightly better than f2 in the senseof the decay of u = eij (t) for large t. For the perturbed version of (3) by the model-implementation anddiﬀerentiation errors, it is also shown that f2 is more robust against those errors than f1 .

4

Conclusion

In this presentation, we formulated rigorously the structure of activation functions in RNN which haveexponential decay of errors and robustness. Moreover, we propose a slightly better activation function
than the super-sigmoid function, which is the best among several candidates discussed in Zhang et al. [2].
Recently, however, Zhang et al. [6] proposed a new approach to find good activation function in RNN.
The argument by taking the new idea into account has not been made, and it is left for further research.

References[1] Y. Zhang and S. S. Ge, ”A general recurrent neural network model for time-varying matrix inversion,”
42nd IEEE Conf. Decis. Cont. pp. 6169-6174, 2003.
[2] Y. Zhang and S. S. Ge, ”Design and analysis of a general recurrent neural network model for timevarying matrix inversion,” IEEE Trans. Neural Netw. vol. 16, no. 6, pp. 1477-1490, 2005.
[3] D. Guo, Z. Nie and L. Yan, ”Novel discrete-time Zhang neural network for time-varying matrix
inversion,” IEEE Trans. Syst., Man., Cybern., Syst., vol. 47 no. 8, pp. 2301-2310, 2017.
[4] Y. Zhang, G. Ruan, K. Li and Y. Yang, ”Robustness analysis of the Zhang neural network for online
time-varying quadratic optimization,” J. Phys. A: Math. Theor., vol. 43, pp. 245-202, 2010.
[5] Y. Zhang, Y. Yang, N. Tan and B. Cai, ”Zhang neural network solving for time-varying full-rankmatrix Moore-Penrose inverse,” Springer Link Comput. vol. 92, no. 2, pp. 97-121, 2011.
[6] Y. Zhang, L. Jin and S. Li, ”Integration-enhanced Zhang neural network for real-time-varying matrix
inversion in the presence of various kinds oh noises,” IEEE Trans. Neural Netw. Learn. Syst. vol. 27,
no. 12, pp. 2615-2627, 2016.

On mode and antimode preserving circular distributions
Toshihiro Abe1 , Yoichi Miyata2 and Takayuki Shiohama3
1
Nanzan University, Japan.
abetosh@ss.nanzan-u.ac.jp
2
Takasaki City University of Economics, Japan.
ymiyatagbt@tcue.ac.jp
3Tokyo University of Science, Japan.
shiohama@rs.tus.ac.jp
Azzalini & Capitanio (2003) constructed an extension of the skew-normal distribution given
by Azzalini (1985) as follows: Let f0 (x) (x ∈ R) be a symmetric pdf about x = 0. Let G be the
distribution function such that G(−x) = 1 − G(x). Then their perturbation of the symmetric
pdf is given by
2f0 (x)G{w(x)},(1)

where the weighting function w : R → R is taken to be an odd function, i.e., w(−x) = −w(x)for x ∈ R. It is noted that there exists an odd function h(x) (−1 ≤ h(x) ≤ 1) such that
2G(w(x)) = 1 + h(x). Then, it is easily shown that∫
(1 + h(x))f0 (x)dx = 1.
RThis point of view is seen in some literature (cf. Huang & Chen, 2007).
Their skew family of the symmetric distributions is not always unimodal even if the basesymmetric distribution f0 (x) is unimodal. On the other hand, Jones (2014) proposed anothermember of skew distributions with pdf
f (x) = f0 (s−1 (x)),

x ∈ Sf ,(2)

where s′ (x) + s′ (−x) = 2, the function s : D → Sf (D ⊇ Sf0 ∋ 0) is a monotone and bijection
function and Sf is a support of the function s. A noteworthy point is that the resulting density
f (x) is always unimodal regardless of the parameter values of the base symmetric unimodal
density.
For a construction of the circular version of the skew distributions by Azzalini & Capitanio(2003), Umbach & Jammalamadaka (2009) adapted the pdf (1) in order to obtain a generalmeans of skewing symmetric circular models. Let f0 (θ) and g(θ) (θ ∈ [−π, π)) be circular pdf’s∫θsymmetric about θ = 0 and G(θ) = −π g(ϕ)dϕ is the distribution function of the latter. The
weighting function w is an odd periodic function with |w(θ)| ≤ π. Thenf (θ) = 2f0 (θ)G{w(θ)},−π ≤ θ < π(3)

is a circular pdf. As a tractable and simple case of the circular pdf (3), Abe & Pewsey (2011)
also introduced the k-sine-skewed circular distributions (SSCDk ) with the densityf (θ) = (1 + λ sin kθ)f0 (θ)
1− π ≤ θ < π,

where λ (−1 ≤ λ ≤ 1) is the skewing parameter. For a related topic for SSCDk with k = 1, Abe& Pewsey (2011) investigated its properties and proposed the sine-skewed Jones-Pewsey family
as well as its three special cases. As they showed, if f0 (θ) is the wrapped Cauchy density, theresulting density is always unimodal. Otherwise, the resulting density is not always unimodal.
Combining with (2) and (3), we obtain another family of unimodal skew distributions on thecircle given byf (θ) = f0 (s−1 (θ; w)),θ ∈ [−π, π),(4)which is an extension of the model in Jones & Pewsey (2012), where f0 is a symmetric unimodal
density about θ = 0 and the function s is given by∫ θs(θ; w) = 2−πG(w(t))dt − π.

(5)

It is not easy to get the explicit expressions for the moments of the family (4) in general,
except for the Jones & Pewsey (2012)’s inverse Batschelet distributions (see Abe, 2015). The
advantage of the family of the distributions is always unimodal. Note that this function is abijection function from [−π, π) to [−π, π), and can be rewritten as s(θ) = θ + H(θ), where H(θ)
is an even function with |H ′ (θ)| ≤ 1. Further, let us confine our attention to the function s
which maps zero to zero. Then the distribution with pdf (4) has the mode-preserving property(cf. Fujisawa & Abe, 2015).

References[1] Abe, T. & Pewsey, A. (2011). Sine-skewed circular distributions. Statistical Papers, 52,
683–707.
[2] Azzalini, A. (1985). A class of distributions which includes the normal ones. Scandinavian
Journal of Statistics, 12, 171–178.
[3] Azzalini, A. & Capitanio, A. (2003). Distributions generated by perturbation of symmetry
with emphasis on a multivariate skew t-distribution. Journal of the Royal Statistical Society,
Series B, 65, 367–389.
[4] Fujisawa, H. & Abe, T. (2015). A family of skew distributions with mode-invariance throughtransformation of scale. Statistical Methodology, 25, 89–98.
[5] Jones, M.C. & Pewsey, A. (2012). Inverse Batschelet distributions for circular data. Biometrics, 68, 183–193.
[6] Umbach, D. & Jammalamadaka, S. R. (2009). Building asymmetry into circular distributions. Statistics & Probability Letters, 79, 659–663.

2

シンポジウム報告書（タイトル:On estimating finite mixtures ofskew-rotationally-symmetric distributions)
宮田 庸一∗ , 塩濱 敬之† , 阿部 俊弘‡

1

はじめに
円周データ, 方向データにおいては, しばしば分布の形が多峰型になったり, モードに関して非対称となることがある. そのようなデータに対しては, 通常回転に対して対称 (rotationally symmetric, 以後単に対称という) な分布の有限混合モデルを用いることが多く, 様々な研究者 (例えば Holzmann, Munk, and Stratmann [3],Banerjee Dhillon, Ghosh and Sra [1]) によりアルゴリスムが提案され, 理論的な性質が明らかにされている. しかしながら, 一方でコンポーネントが非対称な分布の有限混合モデルにおいては, あまり研究がなされていな
い. このため, 本報告においては, 以下の sine-skewed wrapped Cauchy 分布の有限混合モデルを考え, その最尤推定量を求めるための EM アルゴリズムおよび漸近的な性質に関して報告を行った:

f (θ|γ) =g∑

αj fswc (θ|µj , ρj , λj ),

(1)

j=1

ただし g はコンポーネントの数を表し, αj (j = 1, ..., g) は混合比率,

fswc (θ|µ, ρ, λ) =1 − ρ22π(1 + ρ2 − 2ρ cos(θ − µ)){1 + λ sin(θ − µ)}(2)は sine-skewed wrapped Cauchy 分布の確率密度関数を表す. また γj = (µj , ρj , λj )T , γ = (α1 , ..., αg , γ1T , ..., γgT )Tをパラメーターベクトルとする.

2

EM アルゴリズム
大きさ n の観測値 Θ1 ,....,Θn が互いに独立に, sine-skewed wrapped Cauchy 分布の有限混合モデル (1) に従うとする. 我々の提案する EM アルゴリズムにおける E-step は，通常の E-step と同じである. 一方で M-step∗ 高崎経済大学経済学部. E-mail:ymiyatagbt@tcue.ac.jp† 東京理科大学工学部‡ 南山大学理工学部

1

は通常の方法とは異なり, Qj 関数と呼ばれる j 番目のコンポーネントのパラメーターの関数を最大化するために, それぞれのコンポーネントにおけるパラメーター空間の領域を 2 つに分割する. そして, それぞれの領域に(k)

おいて Qj (γj , γ (k−1) ) を最大化する γj の値を求め, 大きい Qj を持つ γj の値を γjることを報告した.

として更新する方法であ

この修正した EM アルゴリスムと通常の EM アルゴリスムをシミュレーションにより比較したところ, 20 回中 8 回はアルゴリスムの収束において改善が見られ, 12 回は通常の方法と同じパフォーマンスであることを確認した. またこのアルゴリスムを実データ (風向データ) に対して適用した例を報告した.

3

一致性
有限混合モデルにおいては, 通常の識別可能性が成りたたないため, 真のモデルを表すパラメーターの集合と, 最尤推定量をプラグインした密度関数の表すパラメーターの集合との”集合間の距離”において一致性を示す必要がある (詳しくは例えば, Cheng and Liu [2] を参照のこと). 本報告においては, この一致性を示すために必要となる, sine-skewed wrapped Cauchy 分布 (2) の識別可能性が成り立つことを報告した. さらにこの結果と Cheng and Liu [2] の定理を用いることで, 集合間の距離の下での最尤推定量の強一致性が成り立つことを
報告した.

4

質疑応答について会場に寄せられたコメントとしては, 我々の提案する EM アルゴリスムの計算時間に関する質問および EM

アルゴリスムの収束に関する質問がされた.

参考文献[1] Banerjee, A., Dhillon, I. S., Ghosh, J. P. and Sra, S. (2005). Clustering on the unit hypersphere using
von Mises-Fisher distributions, Journal of Machine Learning Research, 6, Sep, 1345–1382.
[2] Cheng, R.C.H. and Liu, W.B. (2001). The consistency of estimators in finite mixture models, Scandinavian Journal of Statistics, 28, 4, 603–616.
[3] Holzmann, H., Munk, A. and Stratmann, B. (2004). Identifiability of finite mixtures-with applications
to circular distributions, Sankhyā, 440–449.

2

BOOTSTRAP INFERENCE UNDER RANDOM
DISTRIBUTIONAL LIMITS
Speaker: Giuseppe Cavaliere, University of Bologna(this is a joint work with Iliyan Georgiev, University of Bologna )Asymptotic bootstrap validity is usually understood and established as consistency of the distribution of a bootstrap statistic, conditional on the data,
for the unconditional limit distribution of a statistic of interest. In many applications, however, the bootstrap statistic may possess, conditionally on thedata, a random limit distribution. Among others, cases of random bootstrap
limit distributions are documented for in…nite variance processes (Athreya, 1987;
Knight, 1989; Aue et al., 2008; Cavaliere et al., 2016), time series with unit roots(Basawa et al., 1991; Cavaliere et al., 2015), and parameters on the boundaryof the parameter space (Andrews, 2000). In most cases, the occurrence of arandom limit distribution for a bootstrap statistic given the data –in contrast
to a necessarily non-random limit of the unconditional distribution of the corresponding statistic, computed on the original sample –is taken as evidence offailure of the bootstrap.
In this paper we show that randomness in the limiting distribution of abootstrap statistic need not invalidate bootstrap inference, as the bootstrap
may still deliver con…dence intervals (or hypothesis tests) with the desired coverage probability (or size) when the sample size diverges. Moreover, in suchcases the bootstrap may also have the appealing asymptotic interpretation ofa conditional inferential procedure, and may deliver e¢ ciency (or power) gains
over unconditional inference.
To see why, it is useful to recall that, apart from an unconditional limit distribution, a statistic in general possesses a family of (random) conditional limit
distributions, depending on the choice of the conditioning -algebra. If oneof these conditional (random) limit distributions matches the (random) limit
distributions of a chosen bootstrap statistic, then –under regularity conditions
that will be discussed in the paper –inference based on the bootstrap is asymptotically valid and, importantly, conditional in nature. This observation was
initially made by Lepage and Podgorski (1996) but has not been pursued further in the bootstrap literature, in particular because its development requires
probabilistic tools that are not widely popular in this …eld.
Conditional inference can be justi…ed by the ‘conditionality principle’, according to which “the evidential meaning of any outcome of any mixture experiment is the same as that of the corresponding outcome of the corresponding
component experiment, ignoring the overall structure of the mixture experiment” (Birnbaum, 1962, p.271). Whenever for a statistic of interest the bootstrap estimates consistently a component of the limit unconditional distribution
1viewed as a mixture of conditional distributions, the bootstrap can be regarded
as a large-sample implementation of the conditionality principle.
In such cases the bootstrap replicates asymptotically the property of conditional tests and con…dence intervals to have conditionally constant size andcoverage probability, respectively. Regarding tests, this property has been argued to be necessary for test optimality in the special case of conditioning on
a complete su¢ cient statistic; see Lockhart (2012). More generally, gains in
power and precision can be expected to occur when the reference population ise¤ectively restricted to outcomes that share statistically relevant features withthe actual sample. For instance, in the case of con…dence intervals, Lepage andPodgorski (1996, Figure 2) provide numerical evidence of substantial precision
gains in a particular implementation of a permutation bootstrap with a random
limit of the bootstrap statistic, where conditioning is on the order statistics ofregression residuals (ancillary in that context).
Following a practice in the literature (see, e.g., Lockhart, 2012), we recast theconstant conditional size and coverage probability property into the requirementthat bootstrap p-values should be uniformly distributed conditionally, at least
asymptotically. One of our main results is a general su¢ cient condition forthis to be the case. We also provide conditions for the more basic propertyof unconditional asymptotic distributional uniformity of p-values; this property
implies asymptotic control of the frequency of wrong inferences on average overthe conditioning variables but no longer warrants a conditional interpretationof the bootstrap inferential procedure.
When dealing with random limiting distributions, the usual convergence concept employed to establish bootstrap validity, i.e. weak convergence in probability, can only be employed in some special cases. Therefore, in this paper
we discuss asymptotic bootstrap validity also in cases where consistency of thebootstrap distribution for a conditional (null) limit distribution of an original
statistic holds in a sense weaker than the usual weak convergence in probability.
To show the practical relevance of our results, we include an analysis of thewell-known and much applied (but also misunderstood) bootstrap tests of parameter constancy in regression models where the design matrix could be random
but be conditioned upon; see Hall (1991,p.170). In the resampling process forming the bootstrap sample, it appears natural to take the design matrix as …xed,
i.e. it does not vary across the bootstrap repetitions. Accordingly, bootstrap
algorithms with this feature are sometimes labelled as ‘…xed design’, ‘…xed regressor’ or ‘conditional’ bootstrap. Under a set of assumptions proposed byHansen (2000), we argue that the …xed-regressor bootstrap test statistics haverandom limit distributions, thus invalidating previous claims that the bootstrap
is consistent for the unconditional limit distribution of the original parameter
constancy test statistics. Then we provide conditions under which the …xedregressor bootstrap entails conditional asymptotic inference.

2

Testing for Series Correlation and ARCH Eﬀect ofHigh-Dimensional Time Series Data
Speaker: Shiqing Ling ∗Hong Kong University of Science and Technology(this is a joint work with Ruey S. Tsay in University of Chicago
and Yaxing Yang in Xiamen University)
The small n/large p problem has been extensively studied in the last decadefor independent and identically distributed (iid) data. Various methods have beendeveloped for variable selection in a large data set with a sparse covariance matrix,
suchas the lasso-type estimation procedure in Cai et al. (2011) and random matrix theory in Bai et al. (2015); see also Paul and Aue (2014), and the principal
component procedure in Shen and Huang (2007). There have been several papers
that tried to explore the small n/large p problem in the time series setting. From
example, Song and Bickle (2011) and Han and Liu (2016) proposed two diﬀerent
approaches to estimate the parameters of a large vector autoregressive model.
To make our discussion clear, we ﬁrst deﬁne the high-dimensional time series(HDTS) as follows.
Definition 1.1. A stochastic process Xt = (x1t , · · · , xpt )′ is called a high dimensional (HD) time series if Xt → (x1t , · · · , xit , · · · )′ as p → ∞, where t = 0, ±1,±2, · · · . We say that a HDTS {Xt } is weakly stationary if, for any given p,(a).

EXt = µ(p) ,(b).

E(Xt − µ(p) )(Xt+l − µ(p) )′ = Γl , l = 0, ±1, ±2, · · · ,
(p)

(p)

where µ(p) is a p-dimensional constant vector and Γl is a p × p constant matrix
independent of t. Similarly, we say that a HDTS {Xt } is strictly stationary if it isstrictly stationary any given p. As a special case, if Xt is serially independent of(p)a given p, then Xt is a HD white noise. When Γl = 0 for all l(̸= 0) and p, we
say that the HD time series {Xt } is serially uncorrelated. Let Ft be the σ−ﬁeldgenerated by {Xs : s ≤ t} and denote
E[(Xt − µ(p) )(Xt − µ(p) )′ |Ft−1 ] = Ωt−1 .
(p)∗E-mail address: maling@ust.hk

1

(p)

(p)

The HDTS {Xt } is said to have no ARCH eﬀect if, for any given p, Ωt = Γ0 , a
constant matrix; otherwise, it has ARCH eﬀect.
(p)
For ease in notation, we drop the superscript of µ(p) and Γl in the sequal. Keepin mind, however, that the dimensions of µ and Γl are increasing with p.
Using the random matrix theory, Li, Yao, Lam and Yao (2016) proposed a test forserial correlation of a HD time series. Change, Yao and Zhou (2016) used the maximum absolute auto- and cross-correlation to construct a test for the same purpose,but critical values of the test statistic need to be obtained via a bootstrap procedure.
Recently, Tsay (2017) proposed a rank-based test using Spearman’s rank correlationand the extreme value theory. The asymptotic distribution is derived under the null
hypothesis and contemporaneous independence. The moment condition of the data
is not needed and the asymptotic critical values of test statistics have a close-form
solution. This is a nonparametric approach, see more review on this in Section 3.
Testing for the ARCH eﬀect in a time sereis has become an important issue inanalysis of univariate and multivariate time series since Engle (1982) proposed the
ARCH model. The classical methods for detecting ARCH eﬀect are the Lagrange
multiplier test of Engle (1982) and the Portmanteau test of Li and Mak (1986),
which are based on the squares of random variables. Li and Li (2008) considered
two Portmanteau tests by using the absolute random variables. Zhu and Ling (2011)
proposed a sign based test statistic. Ling and Li (1997) extended Li-Mak test for
multivariate ARCH models by using a quadratic form of random vectors. We refer
to Tsay (2014) for more discussions on the subject. To the best of our knowledge,
there is no test statistic available for detecting the ARCH eﬀect in a HD time series.
This paper proposes two Portmanteau tests for detecting serial correlation and
ARCH eﬀect of a HD time series {Xt } allowing the dimension of data p = p(n) → ∞when the sample size n → ∞. We ﬁrst show that the sample autocorrelation functionof the L1 −norm of data is asymptotically normal and a norm-based Portmanteau
test statistic can be constructed. When the cross-sectional variables are s-dependent(i.e., at most s elements are dependent), the test statistic still works well in the casewith p > n. Using a suitable function of data, the norm-based test can be applied
for the heavy-tailed time series. We next show that the sample rank autocorrelation
function, i.e. Spearman’s rank correlation, of the L1 −norm of data is asymptoticallynormal and the norm-based rank Portmanteau test statistic can also be constructed.
Surprisingly, the norm-based rank test is dimension-free, i.e. independent of p,
and without requiring any moment condition of data or the covariance structure
condition required in the literature. Two standardized norm-based tests are further
discussed. Simulation results show that both test statistics have satisfactory size
and are very powerful even for small n and large p.

2

Cancer Gene Analysis by Singh et al. Microarray Data
成蹊大学 名誉教授 新村秀一（しんむら しゅういち）
1. はじめに
Fisher が計算機のない時代、Fisher の仮説を考えることで Fisher の線形判別関数（F-LDF）を定式化し、判別分析の世界を開き多くの分野で成果を上げた。彼は Iris データで検証し、さらに「現実のデータが仮説を満たさない場合、推測学的結論は信用できないことがある」と明確に述べている。筆者は多くの医学データなどの研究を通して 5 つの問題を見つけた。そして最小誤分類数(Minimum Number of Misclassification, MNM)基準による最適
線形判別関数(Revised IP-OLDF, RIP)[2]と「小標本のための 100 重交差検証法（新手法 1）」と「癌の遺伝子デー
タの Matryoshka 特徴抽出法（新手法 2）」を開発して 6 種の実データと 6 種の Microarray データで実証研究し2017 年 3 月の特任教授退任を前に全ての問題を解決し、昨年末に Springer から解説書を出版した[3]。本発表では、応用研究として行った米国の 6 研究グループが公開している Microarray データ(Dataset)を用いて行った「癌の遺伝子解析から遺伝子診断」の重要点を Singh ら[5]の結果で報告する。
2. Microarray データ
2015 年 10 月 26 日に、「多様な分野における統計科学の新展開」で石井ら[1]の発表で、米国の 6 研究グループ
が Microarray データを集め論文を書いていて、そのデータが HP(http://www.bioinf.ucd.ie/people/ian/)[8]から入手できることを知った。2015 年 10 月 28 日にそれをダウンロードし 8 種の LDF で判別した。3 種の OLDF は MNM=0
で、k 以下の遺伝子の判別係数が 0 でなく残り全て 0 になった。即ち自然に k 個以下の遺伝子で癌遺伝子が特定できる。これを Small Matryoshka(SM1)と呼ぶ。この SM1 を全遺伝子から省いて再度判別すると別の SM2 が求まる。
そして Dataset は複数個の SM の排他的和集合の信号と、高次元の MNM≧１の雑音に分かれる事が分かった。30 年
以上癌の遺伝子解析が行われ、有用な結果が得られなかった（問題 5）。これは、統計手法で雑音を含んだ
Dataset の有効な分析ができない事と、SM の排他的和集合という特殊な構造が分からなかったためと考えられる。
そして、12 月 20 日までに「6 種のデータ全てが、表のように SM の排他的和集合である」ことが分かった。
表

6 個の Dataset の新手法 2 の結果（JMP 列まで）と遺伝子診断（SM 列以降）

章 Dataset2

Alon et al.

3

Alon et al.

4Singh et al.

5

Golub et al.

6

Tien et al.

7

Chiaretti et al.

8

Shipp et al.

2 群と患者数
Normal (22) vs. tumour
cancer (40)
Normal (22) vs. tumour
cancer (40)
Normal (50) vs. tumour
prostate (50)JMP

SM Max Ratio Min Ratio >=5%

PCA

BGS130

0.90%

0.00%

0

5(8.0)

64

26.76%

2.35%

63 30.40%

2(1.6)

179

11.67%

0.28%

38 14.35%

8(11.6)69

15.69%

0.00%

13 34.88%

False (36) vs. True (137) 3（3.9）

159

19.13%

0.63%

27

95

38.98%

10.73%

95 51.46%130

30.67%

4.99%

129 31.70%

All (47) vs. AML (25)

B-cell (95) vs. T-cell10(9.8)(33)
Follicular lymphoma (19)29(16.8)
vs. DLBCL (58)

4.50%

24%

3. 癌の遺伝子解析から癌の遺伝子診断へ
今年に入ってから、遺伝子数が 2000 個と少ない Alon らのデータが筆者が BGS（Basic Gene Set）と呼ぶ癌遺1

伝子 130 組の排他的和集合であることが分かった。BGS は、1 個の遺伝子を省くと MNM が 1 以上になる最小の SM である。Datasets は統計手法でよい結果が出なかったが、各 SM は小標本のため統計的にアプローチできるので、FLDF、2 次判別関数、ロジスティック回帰、一元配置の分散分析、クラスター分析、主成分分析で分析を行った。しかし、ロジスティック回帰だけが BGS と SM のすべてがNM=0 であるが、他の手法は 2 群が線形分離可能な兆候を示さなかった。しかし、BGS と SM を RIP で判別した判別
スコアを新しい変数としたデータを作成した。2 群の平均の t 検定は、t 値が正から、ほぼ 0、そして負のものが得
図
PCA の結果られた。クラスター分析と PCA は 2 群がきれいに分かれた。そのうえ、PCA の第 1 固有値は、青島らの[]PCA による第 1 固有値がスパイク状になる研究と同じである。また、図のスコアプロットのように正常症例は第 1 主成分軸上の負に、がん症例が正にきれいに分かれた。また、線形分離可能な判別分析の評価に SV 間の距離 2 を判別スコ
アの範囲に対する比を RatioSV という統計量で評価することにした。表の章列は、[4]の章番号を示す。6 種の数理計画法による LDF は、SV で-1 以下に正常患者、1 以上に癌患者が判別される。しかし判別結果の評価基準がなかった。SV の距離の 2 が判別スコアの範囲の長さの何%であるかを表す RatioSV（=2/判別スコアの範囲*100）を考えた。これが表の Max Ratio と Min Ratio 列である。Alon 他の 130 個の RIP の中で最大は 0.90%で最小は 0.00%（0.005 未満）である。これで新規の患者の検証標本でよい結果が出るとは考えられない。しかし SM では最大が26.76%で、最小が 2.35%である。この値が 5%以上のものが 63 個あり、1 個だけが 2.35%である。この閾値は今後の課題であるが、仮に 5%以上とすれば、63 個の RIP が癌診断に利用できる。Chiaretti の場合、約 40%の大きな窓が開き残りの 60%に正常と異常が散らばっている。また Ward 法と PCA の結果は綺麗に 2 群に分かれる。図は Alon 他の 64 個の判別スコアデータの PCA の結果である。左の固有値を見ると第 1 主成分（Prin1）の固有値が 39.4 とSpike 状に大きい。これは 2 群が大きく離れ、2 群の群内分散が小さいためと考える。右の因子負荷量は、4 象限から 1 象限に布置している。真ん中のスコアプロットは、正常群が Prin1 の負にほぼ直線上に、癌は原点から正の方
向に扇のように布置している。右の癌患者は悪性で原点近くの癌は比較的軽度の癌患者であり、正常の原点に近い
患者は癌になりやすいか否かは、病理学者が患者カルテを調べれば分ると考えている。即ち、Prin1 の主成分軸は63 個の RIP の判別スコアと同じく癌の悪性度の指標に使えることが期待される。もう一つの応用としては、癌患者が治療で正常領域に誤判別されれば 5 年も経過観察なく治癒したと確定できると素人的に考える。主成分軸上のRatioSV は 30.4％で 26.76%よりも約 4%も大きいのは、64 個の合成のためである。
4．まとめ
残念なことに、せっかく求めた癌の悪性度指標が、公開されたデータの変数が真の遺伝子名を他の変数に置き換えているために医師による検証が行えない。真の遺伝子名を知る方法をご存知の方は、教えていただきたい。
REFERENCES
1.

2.
3.
4.
5.
6.
7.
8.

Ishii A, Yata K, Aoshima M (2014) Asymptotic distribution of the largest eigenvalue via geometricrepresentations of high-dimension, low-sample-size data. Sri Lankan J. Appl. Statist., Special
Issue: Modern Statistical Methodologies in the Cutting Edge of Science (ed. Mukhopadhyay, N.):8194.
新村秀一 (2010). 最適線形判別関数. 日科技連出版.
Shinmura, S. (2016). The New Theory of Discriminant Analysis after R Fisher, Springer
Shinmura, S. (2017). From Cancer Gene Analysis to Cancer Gene Diagnosis, Amazon, 2017年6月.
Vapnik, V. (1995). The Nature of Statistical Learning Theory. Springer-Verlag.
新村秀一 (2011). 数理計画法による問題解決法. 日科技連出版.
新村秀一 (2004).JMP活用 統計学とっておき勉強法.講談社．
Jeffery, IB. et al.(2006). Comparison and evaluation of methods for generating differentiallyexpressed gene lists from microarray data, BMC Bioinformatics. 7:359: 1-16.
2

二元分割表における独立性および対称性のモデルと分解
佐藤 駿
合田 新平
三枝 祐輔
生亀 清貴
山本 紘司

(東京理科大学大学院 理工学研究科)(東京理科大学大学院 理工学研究科)(横浜市立大学 医学部)
(東京理科大学 理工学部)(大阪市立大学大学院 医学研究科)

第一部「既知のスコアを用いた倹約な独立モデルの直交分解」
行と列に順序がある r × c 分割表において, 行変数を X1 , 列変数を X2 とし, P (X1 =i, X2 = j) = pij とする (i = 1, . . . , r; j = 1, . . . , c). 行変数と列変数の独立性が成り立たない場合には, それらの行と列の連関構造に関心がある. 行変数と列変数に既知のスコア{ui }, {vj } をそれぞれ割り振ることができるとき (ただし, u1 < · · · < ur , v1 < · · · < vc ),Tomizawa (1992) は倹約な Linear-by-Linear 連関 (P LL) モデルを次のように提案した:

pij = µαui β vj θui vj(i = 1, . . . , r; j = 1, . . . , c).

特に, {ui = i}, {vj = j} のとき, P LL モデルは倹約な一様連関 (P U ) モデルである.
ここで, g1 (i) = ui (i = 1, . . . , r), g2 (j) = vj (j = 1, . . . , c) とおく. Y1 , Y2 を Y1 =g1 (X1 ), Y2 = g2 (X2 ) によって定義し, Y1 と Y2 の Pearson 相関係数を ρ(Y1 , Y2 ) とする.
また, G2 (M ) をモデル M の適合度を検定する尤度比カイ二乗統計量とする. 本報告では, 行と列に既知のスコアを割り振ることができる場合に, 倹約な独立モデルを定義した. さらに, P LL モデルおよび, ρ(Y1 , Y2 ), Kendall の τb (Kendall, 1945), Spearman のρs (Stuart, 1963; Kendall and Gibbons, 1990, p. 8) を用いて, 倹約な独立モデルの分解について考えた.
定理 1. P I モデルが成り立つための必要十分条件は, P LL モデルと ρ(Y1 , Y2 ) = 0 の両方が成り立つことである.
定理 2. P I モデルが成り立つための必要十分条件は, P LL モデルと τb = 0 の両方が成り立つことである.
定理 3. P I モデルが成り立つための必要十分条件は, P LL モデルと ρs = 0 の両方が成り立つことである.
定理 4. G2 (P I) は, G2 (P LL) と G2 (ρ(Y1 , Y2 ) = 0) の和に漸近的に同等である.

第二部「順序付き正方分割表における非制限正規分布型対称モデル」
行と列が順序のある同じ分類からなる r × r 正方分割表を考える. 正方分割表解析においては, 分類間の対称性の解析に関心があり様々なモデルが提案されている. (i, j) セ
ル確率を pij とする (i = 1, . . . , r; j = 1, . . . , r).

Yamamoto, Nakane and Tomizawa (2016) は次のような制限正規分布型対称 (RNDS)

モデルを提案した:
2

2

pij = µαi +j β i+j γ ij

(i = 1, . . . , r; j = 1, . . . , r).

対称性のモデルが成り立たないときには非対称性の構造に関心があり, Tomizawa(1991) は次のような拡張線形対角パラメータ対称 (ELDPS) モデルを提案した:{
δ j−i γ (j−i)(j+i)/2 ψij (i < j),
pij =
ψij
(i ≥ j),ただし ψij = ψji .
本報告では, RNDS モデルの制限を取り除いたモデルを提案し, そのモデルを用いた
RNDS モデルの分解について述べた.
非制限正規分布型対称 (UNDS) モデルを次のように提案した:
2

2

pij = µα1i α2j β1i β2j γ ij

(i = 1, . . . , r; j = 1, . . . , r).

特に α1 = α2 かつ β1 = β2 のとき UNDS モデルは RNDS モデルである. また, ELDPS
モデルの特別な場合であることに注意する.
行と列が順序のある同じ分類からなる r × r 正方分割表において, 行変数を X1 , 列変数を X2 とし, 平均分散一致 (MV) モデルを次のように定義した:

E(X1 ) = E(X2 ),かつ

V ar(X1 ) = V ar(X2 ),ただし,

E(X1 ) =r∑

kpk· ,k=1

V ar(X1 ) =r∑(k − E(X1 ))2 pk· ,k=1

E(X2 ) =r∑

kp·k ,k=1

V ar(X2 ) =r∑(k − E(X2 ))2 p·k ,k=1

次の定理について述べた:
定理 5．RNDS モデルが成り立つための必要十分条件は, UNDS モデルと MV モデルの両方が成り立つことである．

参考文献
Goodman, L. A. (1979). Journal of the American Statistical Association, 74, 537-552.
Kendall, M. G. (1945). Biometrika, 33, 239-251.
Kendall, M. G. and Gibbons, J. D. (1990). Rank Correlation Methods, 5th edition. Edward
Arnold, London.
Stuart, A. (1963). American Statistician, 17, 23-24.
Tomizawa, S. (1991). Metron, 49, 401-409.
Tomizawa, S. (1992). Biometrical Journal, 34, 129-140.
Yamamoto, K., Nakane, H. and Tomizawa, S. (2016). Journal of Statistics: Advances in
Theory and Applications, 15, 71-84.

正方分割表における対称性の尺度について

東京理科大学 工学部
安藤宗司
東京理科大学 理工学部 田畑耕治
東京理科大学 理工学部 富澤貞男

行と列が同じ分類からなる正方分割表データについて考える．表 1 のデータは，臨床試験における 2 つの治療群（エソメプラゾール群とプラセボ群）のベースライン時点と最終評価時点のカテゴリカルアウトカム（修正 LANZA スコア）に関する変化を表しており，正方分割表データである．修正 LANZA スコアは，“ 0 ” が最も良い評価，“ +4” が最も悪い評価となっている. 表 1 のデータにおいて，プラセボ群よりもエソメプラ
ゾール群の方が治療効果が優れているかどうかに関心がある．
表 1: エソメプラゾール群とプラセボ群におけるベースライン時点と最終評価時点での修正 LANZA スコアの変化；Sugano, Kinoshita, Miwa and Takeuchi (2012).
(a) エソメプラゾール群
最終評価
ベースライン時点
時点
0 +1 +2 +3 +4
0
78 9
26
3
1
+1
1
5
6
4
0
+2
9
1
10
3
1
+3
1
0
1
0
0
+4
3
0
1
1
2
計
92 15 44 11
4

計
117
16
24
2
7
166

(b) プラセボ群
最終評価
ベースライン時点
時点
0 +1 +2 +3 +4
0
41 2
19
0
0
+1
8
0
4
0
0
+2
12 4
14
3
0
+3
0
1
1
3
0
+4
29 7
11
6
0
計
90 14 49 12
0

計
62
12
33
5
53
165

正方分割表データの解析においては，分類間の独立性に代わり対称性や非対称性に関するモデルが用いられる．対称性に関するモデルとして，対称モデル (Bowker, 1948) が提案されている．表 1 のデータに対して対称モデルが成り立つとき，ベースライン時点と最終評価時点の修正 LANZA スコアは同等であることを指している．また，表 1 のデー
タが非対称の構造を持つとき，ベースライン時点に対して最終評価時点の修正 LANZA
スコアは改善または悪化していることを指している．
分割表データにおいて，モデルの適合度を評価する際に，適合度検定統計量（例えば，
尤度比カイ二乗統計量）が用いられる．複数の分割表データに対して対称モデルの当てはまりが悪いとき，非対称性の程度を比較することに関心がある．適合度検定統計量を用いて，複数の分割表データに対して非対称性の程度を比較する際に 2 つの問題点が挙げられる．

1. 分割表のサイズ及びサンプルサイズに依存する
2. 非対称の方向性（改善または悪化等）を区別できない臨床試験において，各治療群のサンプルサイズが等しくない場合もあること，治療効果を比較する際に非対称性の構造を区別する必要があるため，適合度検定統計量を用いて治療群間の治療効果を比較をすることはできない．
順序カテゴリ正方分割表において，Tomizawa, Miyamoto and Hatanaka (2001) は対称性からの隔たりを測る尺度 γ を提案した．この尺度 γ は分割表のサイズ及びサンプル
サイズに依存しないが，非対称の方向性（改善または悪化）を区別できないことに注意する．Tahata, Miyazawa and Tomizawa (2010) は，2 種類の非対称性の方向性を区別することができる尺度 ϕ を提案した．この尺度 ϕ は 2 種類の非対称性の方向性を区別できるが，対称モデルからの隔たりを測っていないことに注意する．
本講演では，分割表のサイズ及びサンプルサイズに依存せず，非対称性の程度と 2 種類の非対称性の方向性を同時に分析することができるベクトル尺度 Φ を提案した．
( )γ; 2 × 1 ベクトル.
Φ=ϕ
提案したベクトル尺度 Φ を用いることにより，2 種類の非対称性の方向性を区別しながら，非対称性の程度を測ることができることを示した．また，提案したベクトル尺度は信頼楕円を用いて，複数の分割表データに対する非対称性の程度を視覚的に比較することに有用であることを示した．

ネットワークメタアナリシスによる Comparative Effectiveness Research と高次漸近理論に基づく推測手法
統計数理研究所

野間 久史

メタアナリシスとは，過去に行われた複数の臨床試験のエビデンスを統合し，総合的な治療効果の評価を行うための方法論であり，Evidence-Based Medicine における重要な研究手法として広く普及している．しかしながら，従来のメタアナリシスの方法
は，原則として，対象となる治療の対比較（1 対 1 の比較）の結果を統合するという単純化されたものであり，「有効性についてのエビデンスが確立された複数の治療の選択肢の中から，どの治療を行うのが最善なのか？（有効性，安全性は？費用対効果は？）」などの問いに必ずしも本質的な答えを与えてくれるものではない．折しも，近年，先進諸国における（超）高齢化社会の到来により，医療費・医療資源の節減および効率的な配分のために，比較効果研究（comparative effectiveness research）が世界的に大きな関心を集めており，このような複数の治療を対象とした有効性・有用性を総合的に評価するための方法論に対する要請は飛躍的に高まっている．ネットワーク
メタアナリシス（network meta-analysis）は，このような背景のもと，急速に普及しつつある，従来のメタアナリシスの方法を一般化し，複数の治療の有効性・有用性を比較・評価することを目的とした方法論である．ネットワークメタアナリシスでは，複数の対象となる治療を含む臨床試験の結果を系統的に集め統合し，治療間の間接比較の情報も併せたエビデンスの統合が行われる．これにより，直接比較の行われていな
い治療間の比較も含めて，対象となった治療法すべての有効性・有用性を比較することが可能となる．
しかしながら、ネットワークメタアナリシスは、この数年で急速に普及した新規な方法論でもあり、現在、実践で普及しているスタンダードな方法でも、さまざまな統計学的問題を抱えている可能性がある。その問題のひとつとして、Brockwell and
Gordon (2001), Noma (2011) などで議論されている推測の妥当性の問題がある。メタ
アナリシスでは、異なる情報源から得られるエビデンスの異質性を考慮した上で情報の統合を行うために、変量効果モデル（random effects model）を用いた解析が行われるのが一般的である。ネットワークメタアナリシスでは、対比に基づく（contrast-based）
方法として、以下の多変量変量効果モデルが広く用いられている。
𝒀𝑖 ~𝑀𝑀𝑀(𝜽𝑖 , 𝑺𝑖 ), 𝜽𝑖 ~𝑀𝑀𝑀(𝝁, 𝚺),

(*)

2
𝑠𝑖1
⎛𝜌 𝑠 𝑠𝑺𝑖 = ⎜ 𝑖12 𝑖1 𝑖2⋮
𝜌𝑖1𝑝 𝑠𝑖1 𝑠𝑖𝑖⎝

𝜌𝑖12 𝑠𝑖1 𝑠𝑖22𝑠𝑖2⋮
𝜌𝑖2𝑝 𝑠𝑖2 𝑠𝑖𝑖

⋯
…
⋱
⋯

𝜌𝑖1𝑝 𝑠𝑖1 𝑠𝑖𝑖
𝜏12⎛𝜅 𝜏 𝜏
𝜌𝑖2𝑝 𝑠𝑖2 𝑠𝑖𝑖 ⎞, 𝚺 = ⎜ 12 1 2
⎟
⋮
⋮
2
𝑠𝑖𝑖
𝜅1𝑝 𝜏1𝜏𝑝⎠⎝

𝜅12 𝜏1 𝜏2
𝜏22⋮
𝜅2𝑝 𝜏2 𝜏𝑝⋯ 𝜅1𝑝 𝜏1 𝜏𝑝… 𝜅2𝑝 𝜏2 𝜏𝑝 ⎞
⎟
⋱
⋮
⋯
𝜏𝑝2
⎠

𝒀𝑖 は、i 番目の試験における、特定の参照レベル（プラセボなど）に対しての対象となる治療の治療効果の指標（ハザード比，オッズ比など）の推定量からなる確率ベク
トルであり、𝜽𝑖 はその真値を表すベクトル、𝑺𝑖 は試験内共分散行列（𝒀𝑖 の共分散行列）、

𝚺 は試験間共分散行列を表す。パラメータの推定は、制限付き最尤法（restrictedmaximum likelihood; REML）によって行われるのが一般的である。通常、試験間の異質性を表す 𝚺 に関する情報は、𝒀𝑖 のデータ数に比例することとなり、試験数が十分に大きいもとで、大標本理論による正当化が可能となる。しかしながら、多くのメタ

アナリシスでは、統合の対象となる試験数が十分に大きくはなく、古典的な 2 つの治
療法の対比較のメタアナリシスでは、20 にも満たない試験での統合が行われるのが一般的である。この場合、大標本近似の崩れから、信頼区間の被覆率が名目水準を大きく下回るなどの問題があることが知られている（Noma, 2011）。ネットワークメタア
ナリシスでは、複数の治療法の比較を対象としていることもあり、相対的に対象となる試験数は大きくなるが、サブグループ解析などの感度解析なども含めれば、数試験
程度の条件下で統合解析を行う状況は多く存在する。
本講演では、シミュレーション実験を通して、まず、現状のスタンダードな推測手法である REML 法に、上記のような条件下で、推測の妥当性が成り立たない状況が多く存在することを示す。そして、演者らが近年開発した、ネットワークメタアナリシ
スの解析モデルに汎用的に用いることができる高次漸近理論に基づく有効な推測手法（Noma et al., 2017）の紹介を行う。また、シミュレーション実験および統合失調症などのネットワークメタアナリシスの事例解析を通して、その実践的有用性について示す。
文献
Brockwell SE, Gordon IR. A comparison of statistical methods for meta-analysis. Stat Med.
2001;20(6):825-840.
Noma H. Confidence intervals for a random-effects meta-analysis based on Bartlett-type corrections. Stat
Med. 2011;30(28):3304-3312.
Noma H, Nagashima, K, Maruo K, Gosho, M., Furukawa TA. Bartlett-type corrections and bootstrap
adjustments of likelihood-based inference methods for network meta-analysis. Stat Med. In Press, 2017.