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

多様な高次元モデルにおける理論と方法論，及び，関連分野への応用
基盤研究 (S) 18H05290 (2018 年度〜2022 年度) 「広汎な観測に対する因果性の導入

科学研究費補助金とその最適統計推測論の革新（研究代表者: 谷口 正信）」によるシンポジウムを下記のように催しますので，ご案内申し上げます．

誠(筑波大学)

矢田和善(筑波大学)

石井(東京理科大学)

青嶋晶

記
日時：2020 年 2 月 17 日 (月)〜18 日 (火)場所：イーアスホール(イーアスつくば 2 階 https://tsukuba.iias.jp/floor/floor2.html )〒305-0817 茨城県つくば市研究学園 5 丁目 19 番地

プログラム
2 月 17 日 (月)13:30〜13:40

開会

13:40〜14:20

藤森洸

（早稲田大・基幹理工学部）

高次元・スパースな設定における確率過程の統計モデルに対する Dantzig selector

14:30〜15:10

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

Sangyeol Lee (Department of Statistics, Seoul National University)
Asymptotic properties of mildly explosive processes with locally stationary disturbance
15:25〜16:05

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

多変量バーンバウムサンダース型分布: 非負データの密度推定への応用

16:15〜16:55

本田敏雄（一橋大・経済学研究科）

The de-biased group Lasso estimation for varying coeﬃcient models
17:40〜懇親会2 月 18 日 (火)10:40〜11:10

小西啓介（筑波大・数理物質科学研究科）

矢田

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

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

High-dimensional covariance matrix estimation under the strongly spiked eigenvalue model
11:15〜11:45

江頭健斗（筑波大・数理物質科学研究科）

矢田

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

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

Asymptotic properties of distance weighted discrimination and its bias correction
in HDLSS settings
11:50〜12:20

藤澤健吾（東京理科大・理工学研究科）

田畑耕治（東京理科大・理工学部）

順序カテゴリをもつ多元分割表における非同等性のモデリング12:20〜13:45
13:45〜14:25

昼食
石井

晶

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

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

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

データ変換を用いた高次元次判別分析について

14:35〜15:15

菅澤翔之助（東京大・空間情報科学研究センター）

Robust bayesian modeling with synthetic posterior
15:30〜16:00

桃崎智隆，齋藤 健

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

三枝祐輔（横浜市立大・医学部）

生亀

清貴（日本大・経済学部）

中川智之，石井 晶，富澤貞男

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

分割表における種々の対称モデルに関する尺度

16:05〜16:45

島谷健一郎（統計数理研究所）

荒木希和子（立命館大・生命科学部）

地上部と地下部のデータを統合させた推移行列モデル

16:45〜16:50

閉会高次元・スパースな設定における確率過程の統計モデルに対する Dantzig selector
藤森洸
早稲田大学 基幹理工学部数学科

次の線型回帰モデルを考える．

Y = Zβ + ϵ,
ただし，Y ∈ Rn は目的変数，Z = (Z1⊤ , . . . , Zn⊤ )⊤ ∈ Rn×p は計画行列，ϵ = (ϵ1 , . . . , ϵn )⊤ は平均 0 の独立
同分布な確率変数からなる，Z とは独立な確率ベクトル，β ∈ Rp は未知の回帰係数である．従来の回帰分析
においては，データ (Y, Z) に基づいて，パラメータ β の推定を最小二乗法，あるいは最尤法に基づいて行ってきた．特に，誤差項が正規分布に従うとき，適当な正則条件の下で上記の手法により得られた推定量が漸近的に良い性質を持つことはよく知られた事実である．
一方で，本講演ではパラメータ β ，あるいは説明変数 Zi , i = 1, . . . , n の次元 p がサンプル数 n に比べて大きく，パラメータの非零成分の個数が小さいことを意味する，「高次元・スパース」と呼ばれる設定の下で回帰モデルの推定問題を扱う．高次元の設定は，目的変数に回帰しうる説明変数の成分がサンプル数に比べても多いという状況を想定しており，スパース性の設定は目的変数に真に影響を及ぼす説明変数の個数が少ない，
すなわち説明変数の成分の中に無駄な情報が含まれているという状況を想定している．高次元の設定における最小二乗推定量は計算可能性や，一致性をはじめとする漸近的な良さが保証されないケースが多い．また，ス
パースの設定においてはパラメータの推定量に求められるのはノルムの意味での一致性だけではなく，目的変数に真に回帰する説明変数を正しく選択していることを意味する，変数選択の一致性である．このような背景から，一致性と変数選択の一致性を同時に満たす推定量の構成に関する研究が 1990 年代から盛んに行われてきた．特によく知られているのは Lasso 推定量に代表される，罰則項付き推定関数により構成される推定量である．例えば，線型回帰モデルに対する Lasso 推定量は次の凸最適化問題の解として定義される．

Deﬁnition 1. 次で β の推定量を定義する．
{
β̂L = arg minβ}1∥Y − Zβ∥22 + λ∥β∥1 ,nただし，λ ≥ 0 は解析者が設定する tuning parameter である．
ここで推定に用いている推定関数は，従来用いられてきた最小二乗誤差関数に，パラメータの l1 ノルムで定義された罰則関数を加えたものである．もし，λ = 0 であれば，得られる推定量は最小二乗推定量に一致し，λ を大きく取れば罰則関数 λ∥β∥1 の影響が大きくなって，推定量にスパース性が誘導される．適当な正則条件の下で推定関数を決定する tuning parameter の値を適切に設定することにより，得られる推定量は従来の意味での漸近的な良さと変数選択の一致性を両立する性質を持つことが知られている．さらに，上記の手法は「最小二乗誤差」の部分を一般のモデルにおける従来の M-推定に用いる推定関数に，罰則関数の部分をl1 ノルムそのものではなく，ノルムに関する多様な罰則関数に拡張することができるため，多様な推定手法が提案され現在も盛んに研究されている．
一方，本講演で扱うのは罰則項付き推定問題とは異なるアプローチで構成される，Dantzig selector と呼ばれる推定量である．例えば，線型回帰モデルに対する Dantzig selector は次の不等式制約付き最適化問題の解として定義される：

Deﬁnition 2. 次で β の推定量を定義する．
β̂D := arg min ∥β∥1 ,β∈CC := {β ∈ Rp : ∥ψn (β)∥∞ ≤ λ} ,ただし，ψn は次で定まるスコア関数；

ψn (β) :=1 ⊤Z (Y − Zβ)nであり，λ ≥ 0 は tuning parameter である.

λ = 0 の場合には不等式制約の部分が等式になり，得られる推定量は従来の最小二乗推定量と一致する．また，λ > 0 の場合，不等式制約により最小二乗推定量の近傍が定まり，l1 ノルムによりスパース性が誘導された推定量が得られる．適当な正則条件下で，λ を適切に選ぶことにより，β̂D は Lasso 推定量と同様，漸近的な良い性質と変数選択の一致性を満たすことが知られている．さらにこの手法は，不等式制約を課す関数をZ-推定で用いるスコア関数に置き換えることにより，一般のモデルに対する拡張を考えることは容易である．
また，スコアがパラメータに関して線型関数になるモデルに関しては，推定量を定義する最適化問題が線型計画法に帰着できるということが知られており，実装の意味で利点を持つことが知られている．

Dantzig selector 型の推定量は Lasso 型の推定量とは別のアプローチで構成されたものであるが，得られる推定量の漸近的な性質は互いに似通っている．しかし，Lasso 型の推定量を様々なモデルに拡張する研究が盛んに行われているのに対して，Dantzig selector 型の推定量を拡張する研究は未開拓の部分が多い．そこで本講演では上記の Dantzig selector を確率過程，特に拡散過程の統計モデルに適用し，得られる推定量の漸近的性質に関する議論を行った．なお，本講演の内容は下記の Fujimori (2019) に基づくものである．

参考文献
Fujimori, K. The Dantzig selector for a linear model of diﬀusion processes. Stat. Inference Stoch. Process.
22, no.3, p.475-498. (2019).

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







 

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

#'

¾

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

½+      

 





½





½ '#

½

¾

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





#

 

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

            ,                ,(  !( 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









Æ-多変量バーンバウムサンダース型分布: 非負データの密度推定への応用
柿沢 佳秀 (北大経済)1. はじめに カーネル密度推定量の漸近的性質はよく知られている (Silverman(1986),
Wand and Jones(1995)) が, 推定される密度 f の台が Rd でない (典型的に [0, 1]d や[0, ∞)d の台) のときに, 境界付近で O(1) のバイアスがあり, 一致性をもたないという境界バイアス問題がある. 事実, カーネル密度推定の背後にある畳み込み近似




1
x −v
f (v) dv = Rd k(u)f (x − hu) du ≈ f (x) as h → 0 (カーネル k に
Rd h d khk(−v) = k(v), Rd k(v) dv = 1 を仮定する) が, 有界閉領域や正象限の境界付近では成立しない. この境界バイアスの対処として, 台の変換法 (Marron and Ruppert(1994))の他, renormalization, reflection, 一般化ジャックナイフのような境界カーネルへと修正するアプローチ (Jones(1993)) が主に１次元データで議論されてきた. 報告者はChen(1999,2000) 以降の『非対称カーネル法 (AK 推定法)』に焦点を絞る.
1.1. 1 次元の AK 推定 1 次元の非負データ X1 , . . . , Xn ∼ iid f に基づいた (点x ∈ [0, ∞) = R+ での)AK 密度推定量は n−1 ni=1 kx,b(Xi ) と定義される1 . ここに,
b = b(n) > 0 は n → ∞ のとき適当な速さで 0 に収束する平滑化パラメータ, kx,b(·) は
f と同じ台 R+ を持つカーネルである2.
1.2. 多次元の (積型の)AK 推定 d 次元の非負データ X1 , . . . , Xn ∼ iid f (ただし,Xi = (Xi1 , . . . , Xid ) とする) に基づき f (x), x = (x1 , . . . , xd ) ∈ Rd+ を推定するなら,
いわゆる “積カーネル法”が簡単で, カーネル kx,b (·) を 1.1 節のように選び, 積カーネル
を作り, 積型の AK 密度推定量 n−1 ni=1 dj=1 kxj ,b (Xij ) を定義する. Bouezmarni andRombouts(2010) がベータ/ガンマカーネルによる積型の AK 推定を考察し, 他の文献もある (Funke and Kawka(2015), Zougab et al.(2018), Igarashi and Kakizawa(2020)).
2. BS 型分布 本報告では, 非積型の AK 推定を開発するために, BS 分布の周辺から分布論を再検討した.


A(s/β), s ∈ R+
2.1. 準備 (分布論) １次元の確率変数 ξ が密度関数 p(s; α, β) = φ a(s/β)ααβ(Birnbaum and Saunders(1969)) を持つとき, ξ は BS 分布 BS(α,β) に従う. ここに


,
1
1
1
1
u2
1/2
√
α, β > 0 であり, a(u) = u − u1/2 , A(u) = (1/2) u1/2 + u3/2 , φ(u) = 2π exp − 2 .
また, [a] 非心拡張, [b] 正規分布ベースから対称分布ベースへの拡張, [c] BS と LN とを繋ぐような拡張 (ジョンソンシステムの 1 例とみなせる), [d] 歪分布ベースへの拡張もできて, BS 型分布の一般形は以下のように与えられる.
q-BS 型密度 (柿沢 (2017)): h を R 上の分布として,q aq (s/β)p(h) (s; α, β, θ) = h
1

α

−θ

 Aq (s/β)αβ(s ∈ R+ )この形式は “重み付き推定量” (例えば, Siliverman(1986, Subsection 2.9: General weight function
estimators)) であるが, Chen(1999,2000) 以降は “非対称カーネル推定量”がキーワードである.
2
先行研究 (Jin and Kawczak(2003), Scaillet(2004), Koul and Song(2013), Marchant et al.(2013),
Saulo et al.(2013) など) のオリジナルの推定量において(i) f(0) = 0 の制約 (明らかに, f (0) > 0 の場合には適さない)(ii) O(1) の境界バイアス (Rosenblatt–Parzen カーネル法と同様, これを回避する提案をすべきである)(iii) 積分分散の主要項の存在性に関わる問題 (その存在を仮定してもよいが, f (0) > 0 なる [0, ∞) 上の有界な密度関数を排除してしまう; この仮定を外すべきだろう)などがあり, いずれかの問題があるものについては, それらが回避された “再定義”を用いるべきである.

(この q-BS 型密度は, q = 0 として対数対称密度や対数 skew 密度を含む). ここに,
⎧
⎪
⎨

1 q
(u − u−q ) , q = 0 ,2q
a (u) =,⎪⎩ log u ,q=0
qAq (u) = (d/du)aq (u) .

[例] h(u) = φ(u) の場合が, 正規分布ベースの q-BS, fSg (u) = Cg g(u2) の場合が,対称分布ベースの q-BS である (なお, Dı́az-Garcı́a and Leiva(2005) は θ = 0 の中心ケースで q = 1/2 を扱い, “GBS”と呼んだが, Kakizawa(2018) は一般化の意味を強調した別名から “対称分布ベースの q-BS”を採用). 柿沢 (2017) は, h として Azzalini型uskew/epsilon-skew 密度 hASg,ġ,λ (u) = 2fSg (u)FSġ (λu), hESg, (u) = fSg 1−sgn(u) (u ∈ R)
t
を 2 例としている. ここに, FSg (t) = −∞fSg (u) du (t ∈ R), 及び, sgn(u) = u/|u| if
u = 0 and sgn(0) = 0 とする.
2.2. 多次元拡張 (2 次元の場合のみを掲載) h を R2 上の分布として, aq1  (s /β )qp(h) (s; α, β, θ) = h

1

1

α1

− θ1 ,

 2 Aqj  (s /β )aq2  (s2 /β2 )
j
j
− θ2
α2
α
β
j j
j=1

(s ∈ R2+ ) .

本報告では, h として楕円密度

hEg,ρ (u1 , u2 ) =u21 + u22 − 2ρu1 u2 Cg,2g(u ∈ R2 )2
1/2
2(1 − ρ )2(1 − ρ )を例とし, それをベースにした (q1 , q2 )-BS 分布を扱った.

3. 密度推定への応用 想定したカーネル (報告で, 楕円分布ベースの q-BS 型に限定)に付随するパラメータを推定位置に依存させ, ノンパラメトリックな (可変的)AK 密度


1/2

b
推定量を構築できる. パラメータ (α, β, θ) に対し, αb,c (x) = x+bc
, βb,c (x) = x + bc,
θb,c (x) = θ0 αb,c (x) (ただし, c > 0, θ0 ∈ R) を採用し, １次元の対称分布ベースの q-BSq
カーネル密度推定量 fb,g (x) = n1 ni=1 pgq (Xi ; αb,c (x), βb,c (x), θ0 αb,c (x)), x ∈ R+ (中心ケース θ0 = 0 は Kakizawa(2018), 非心ケース θ0 = 0 は柿沢 (2017) の Azzalini 型
skew/epsilon-skew q-BS カーネルで “歪みパラメータ”を 0 にした特別な場合) を定義したのと同様にして, 積型/非積型3の楕円分布ベースの (q1 , q2 )-BS カーネル密度推定量(d 次元の推定量へ拡張可能)



n 2
1
pqj  (Xij ; αb,c (xj ), βb,c (xj ), θ0 αb,c (xj )) ,n i=1 j=1 g
x = (x1 , x2 ) ∈ R2+ (1)n1q1 ,q2 q ,q fb,g,ρ
(x) =p 1 2 (Xi ; αb,c (x), βb,c (x), ρ, θ0 αb,c (x)) ,n i=1 g,2

q1 ,q2 

fb,g

(x) =

を提案した. ここに, αb,c (x) = (αb,c (x1 ), αb,c (x2 )) , β b,c (x) = (βb,c (x1 ), βb,c (x2 )) とする.
主結果として, (1) の漸近バイアス/分散, MSE,MISE, 強一致性, 漸近正規性を示した(いくつかの主要項は q = (q1 , q2 ) に依存せず, これは初期の正規分布ベースの q = 1/2の BS カーネル密度推定量と q = 0 の LN カーネル密度推定量が, 多くの同じ漸近性質を持つ根拠であった).
3

楕円分布においては, 多変量正規分布以外で, 相関行列を仮に単位行列(無相関) としても, 各要素は

q1 ,q2 
q
,q

1
2
独立にならないことが知られており, 従って, 一般には, f= fとなる. そのため, “積型のb,g

AK 密度推定量”も個別に扱う.

b,g,0

The de-biased group Lasso estimation for varying
coeﬃcient models
一橋大学大学院経済学研究科

本田敏雄

近年データ収集技術の進歩により，様々な分野において非常に多くの高次元デー
タが利用可能になり，それらのデータ解析の必要性が高まってきた．Lasso は Tib-shrani(1996) で提案された，高次元線形回帰モデルに対する推定法であり，高次元
データ解析の必要性，重要性の高まりとともに，実用面および理論面での研究が進んできた．理論面では，Bickel et al.(2009) が最も重要な論文である．
研究の進展に伴い Lasso の問題点：１．特殊な条件下を除けばモデル選択に関する一致性はない，２．比較的バイアスの大きな推定値を与える，が明らかになってきた．後者の問題を解決するために Javanmard and Montanari(2014)， van de Geer etal.(2014)，Zhang and Zhang(2014) により提案されたのが de-biased Lasso である．厳
しい仮定の下ではあるが，Lasso のバイアスを取り除き，高次元の設定を維持した
ままでの統計的推測を可能にした．

Honda(2020) では，de-biased Lasso の手法を高次元変動係数モデルに適用し，理論的な結果およびシミュレーションの結果等を与えている．本講演では Honda(2020)に基づき，変動係数モデルの場合の de-biased Lasso 推定の理論的な結果とシミュレー
ションの結果を報告した．
まず以下のように de-biased Lasso の理論を簡単に紹介した．
高次元線形モデル : X = (Xi j ) を n × p のデザイン行列とし， p > n とする．

Y = Xβ0 + ϵ, ここで Y ∈ Rn および β0 ∈ R p .
Lasso 推定量 : 上のモデルでは，0 でない係数は p に比べて非常に少ないとする．
Lasso 推定量は以下で定義される．
p
∑
{
}
−1
2
b
β B argmin n ∥Y − Xβ∥ + 2λ|β j | ,
β∈R p

j=1

ここで β = (β1 , . . . , β p )T とし， S B { j | β j , 0} かつ s B |S| である．

1

b は以下の式を満たす．
この βb + λκ,n−1 X T Y = n−1 X T X β
ここで κ = (κ1 , . . . , κ p )T ∈ R p および |κ j | ≤ 1 である．
上記の設定では， bΣ = n−1 X T X の逆行列が存在しないので，その代わりとなる

b を Lasso などを用いて構成する．そして以下の式により，de-biased Lasso 推定量
Θは定義される．

bb + Θλκb + 1 ΘXbb =βb T (Y − X β)
bBβ
n
1b T
= β0 + ΘX
ϵ + ∆,n
ここで ∆ は無視可能な項である．
ついで変動線形モデルを紹介した上で，以上の方法に倣って変動線形モデルに対する de-biased Lasso 推定量を定義し，Honda(2020) で与えられている結果を紹介
した．
参考文献

P. J. Bickel, Y. Ritov, and A. B. Tsybakov (2009). Simultaneous analysis of Lasso and
Dantzig selector. Ann. Statist., 37:1705–1732.
A. Javanmard and A. Montanari (2014). Confidence intervals and hypothesis testing forhigh-dimensional regression. J. Machine Learning Research, 15:2869–2909.
T. Honda (2020). The de-biased group Lasso estimation for varying coeﬃcient models.
Forthcoming in Ann. Inst. Statist. Math.
R. Tibshirani (1996). Regression shrinkage and selection via the lasso. J. Royal Statist.
Soc. B, 58:267–288.
S. van de Geer, P. Bühlmann, Y. Ritov, and R. Dezeure (2014). On asymptotically optimal
confidence regions and tests for high-dimensional models. Ann. Statist., 42:1166–1202.
C.-H. Zhang and S. S. Zhang (2014). Confidence intervals for low dimensional parameters in high dimensional linear models. J. Royal Statist. Soc. B, 76:217–242.

2

High-dimensional covariance matrix estimation
under the strongly spiked eigenvalue model
Keisuke Konishia , Kazuyoshi Yatab , Makoto Aoshimaba Graduate School of Pure and Applied Sciences, University of Tsukubab Institute of Mathematics, University of Tsukuba

1

Introduction

One of the features of high-dimensional data is that the data dimension d is high, however, thesample size n is low. This is the so-called “HDLSS” or “large d, small n” data. For HDLSS
data, the sample covariance matrix does not have its inverse matrix. The estimation for the
inverse matrix of a covariance matrix is a crucial issue for high-dimensional data analyses,
especially for pathway analysis and graphical modeling.
Bickel and Levina (2008, AS) gave a thresholding estimator for the inverse matrix of acovariance matrix when the covariance matrix is sparse and its eigenvalues are bounded. However, such sparsity conditions are severe for actual data and often out of touch with reality.
In fact, Aoshima and Yata (2018, Sinica; 2019, MCAP) showed that the first several eigenvalues diverge as d grows and the bounded-eigenvalues condition is quite strict for microarraydata sets. Fan et al. (2013, JRSS-B) proposed a diﬀerent thresholding estimator called theprincipal orthogonal complement thresholding (POET) under the assumption that the first
several eigenvalues diverge rapidly at the rate of d. Unfortunately, the assumption requiredin the POET cannot express the structure of actual data.
Suppose we have a d × n data matrix X = [x1 , . . . , xn ], where xi , i = 1, . . . , n (< d),are independent and identically distributed (i.i.d.) as a d-dimensional distribution with mean
zero and covariance matrix Σ. We denote the eigen-decomposition of Σ by Σ = HΛH T ,where Λ = diag(λ1 , . . . , λd ) with eigenvalues λ1 ≥ · · · ≥ λd ≥ 0 and H = [h1 , . . . , hd ]is an orthogonal matrix with eigenvectors h1 , . . . , hd corresponding to the λ1 , . . . , λd . Let
xj = HΛ1/2 z j , where z j = (z1j , . . . , zdj )T is considered as a sphered data vector having thezero mean vector and identity covariance matrix. We assume that4 ) < ∞ for all r, and(C-i) lim supd→∞ E(zri2 z 2 ) = E(z 2 )E(z 2 ) = 1 and E(z z z z ) = 0 for all r ̸= s, t, u.
E(zriri si ti ui
si
ri
si

When xj s are Gaussian, (C-i) naturally holds.
The sample covariance matrix is given by S = n−1 XX T . Let λ̂1 ≥ · · · ≥ λ̂d ≥ 0 be the∑Teigenvalues of S. Then, we denote the eigen-decomposition of S by S = di=1 λ̂i ĥi ĥi , whereĥi is a unit eigenvector corresponding to the λ̂i . The dual sample covariance∑n matrixT is given
T−1
by S D = n X X. We have the eigen-decomposition of S D by S D = i=1 λ̂i ûi ûi , whereûi is a unit eigenvector corresponding to the λ̂i . We assume the following spiked model forthe eigenvalues of Σ:1

λi→ ∞ as d → ∞ for i = 1, . . . , m, and λi ∈ (0, ∞) as d → ∞ for all i ≥ m + 1.
d1/2
Here, m is a positive and fixed integer.

(C-ii)

Note that (C-ii) is one of the strongly spiked
models given∑ eigenvalue∑dby AoshimaTand Yata
T(2018, Sinica). We devide Σ into Σ1 = m
λ
h
h
and
Σ
=
2
j=m+1 λj hj hj , so thatj=1 j j j
Σ = Σ1 + Σ2 .

2

NR-POET and its asymptotic properties

In this section, we propose a new estimation of Σ−1 . Yata and Aoshima (2012, JMVA)
proposed an eigenvalue estimation called the noise-reduction (NR) method as followings.
∑
tr(S D ) − jl=1 λ̂l
λ̃j = λ̂j −, h̃j = (nλ̃j )−1/2 X ûj (j = 1, ..., n − 1).
(1)n−j
For estimating xTl hj , Aoshima and Yata (2018, Sinica) showed that xTl ĥj and even xTl h̃jinvolve a huge bias. In order to overcome the inconvenience, they gave the modofied NR
method. We modify xTl h̃j as xTl h̃jl byh̃jl =( n ) X ûn1/2 X ûjl
jl
=
,
1/2
n − 1 (nλ̃j )1/2
(n − 1)λ̃whereûjl = (ûj1 , ..., ûjl−1 , 0, ûjl+1 , ..., ûjn )T .

j

We propose a new estimation of Σ−1 by applying the noise-reduction (NR) method to the∑T
POET. By using the modified NR method, we estimate wl by w̃l = (I d − m
j=1 ĥj h̃jl )xl .
( )e 2 , where Σ
e 2 = ∑n w̃l w̃T /n and T (·) is sparse operator
We consider estimating Σ2 by T Σ
l
l=1
( )
e 2 byproposed by Bickel and Levina (2008, AS). We denote the eigen-decomposition of T Σ( ) ∑d
( )
T
e2 =
e
T Σ
j=1 λ́j h́j h́j , where λ́j s are eigenvalues of T Σ2 having λ́1 ≥ · · · ≥ λ́d ≥ 0 andh́j is a unit eigenvector corresponding to the λ́j . Note that
Σ−1 =

m
∑

T
λ−1
j hj hj +

j=1

d
∑

T
λ−1
j hj hj .

j=m+1

Finally, by applying the NR method to the intrinsic part and the noise part, we propose to
estimate Σ−1 by
(
)(d−m
)(
)
m
m
m
∑
∑
∑
∑
−1
T
T
T
T
−1
−1
e
Σ =
λ̃j h̃j h̃j + I d −h̃j h̃jλ́j h́j h́jId −h̃j h̃j .
j=1

j=1

j=1

j=1

Note that the intrinsic and noise parts are orthogonal. We call this new estimation methodthe ”NR-POET”. We have the following result.
Theorem 2.1. Assume (C-i) and (C-ii). Then, under some regularity conditions, it holds
that as d → ∞ and n → ∞√e −1 Σ1/2 − I d ||F / d = oP (1).
||Σ1/2 Σ
2Asymptotic Properties of Distance Weighted
Discrimination and Its Bias Correction
in HDLSS Settings
Kento Egashiraa , Kazuyoshi Yatab , Makoto Aoshimaba1

Graduate School of Pure and Applied Sciences, University of TsukubabInstitute of Mathematics, University of Tsukuba

Introduction

Suppose we have two independent and d-variate populations, Πi , i = 1, 2, having an
unknown mean vector µi and unknown covariance matrix Σi (> O) for each i = 1, 2. Let
us have a trainning data set, xi1 , . . . , xini , from Πi for i = 1, 2. Let N = n1 +n2 . We denotethe class labels of yi by −1 for i = 1, . . . , n1 , and by +1 for i = n1 + 1, . . . , N . We simply
write that (x1 , . . . , xN ) = (x11 , . . . , x1n1 , x21 , . . . , x2n2 ). Let x0 be an observation vectorof an individual belonging to one of the Πi s. We assume that x0 and xij s are independent.
The support vector machine (SVM) developed by Vapnick (2000) is well known as arepresentative tool of binary linear discriminant analysis. However, Marron et al. (2007)
pointed out that the SVM causes data piling in the HDLSS context. Data piling is aphenomenon that the projection of a trainning data to the normal direction vector of aseparating hyperplane is same for each class. In order to avoid the data piling problemof the SVM, Marron et al. (2007) proposed distance weighted discrimination (DWD).
Whereas the SVM ﬁnds the optimal hyperplane by maximizing the minimum distancesfrom each class to the hyperplane, the DWD ﬁnds a proper hyperplane by minimizing thesum of reciprocals of the distance from each data point to the hyperplane. The DWD
cares all the data vectors that are not always used in the SVM. Unfortunately, the DWD
is designed for balanced trainning data sets. For imbalanced trainning data sets, Qiao
et al. (2010) developed weighted DWD (WDWD) that imposes diﬀerent weights on twoclasses. However, the WDWD is sensitive for a choice of weight for each class and it isquite diﬃcult to select a proper weight in actual data analyses.
In this talk, we investigated DWD theoretically in the HDLSS context where d → ∞while n is ﬁxed. We showed that the DWD includes a huge bias caused by heterogeneityof covariance matrices as well as sample imbalance. We proposed a bias corrected-DWD(BC-DWD) and showed that the BC-DWD can enjoy consistency properties about misclassiﬁcation rates.

1

2

DWD and its asymptotic properties

We deﬁne the classiﬁer for the DWD as ŷ(x0 ) = ŵT x0 + b̂. One classiﬁes x0 into Π1
if ŷ(x0 ) < 0 and into Π2 otherwise. Let ∆ = ∥µ1 − µ2 ∥2 . We assume the following
assumptions:Var(∥xij − µi ∥2 )(A-i)
→ 0 as d → ∞ for i = 1, 2;∆2
tr(Σ2i )→ 0 as d → ∞ for i = 1, 2.
∆2Let ∆∗ = ∆ + tr(Σ1 )/n1 + tr(Σ2 )/n2 and(A-ii)

3/2

δ=−

3/2

(n1 − n2 )∆∗tr(Σ1 ) − tr(Σ2 )+.
√√N ( n1 + n2 )∆N∆

We consider the following assumption:mini=1,2 ni(A-iii) lim sup |δ| <.
Nd→∞Let e(i) denote the error rate of misclassifying an individual from Πi into the other classfor i = 1, 2.
Theorem 1. Under (A-i) to (A-iii), the DWD holds:e(1) → 0 and e(2) → 0 as d → ∞.
Corollary 1. Under (A-i) and (A-ii), the DWD holds:
n2
; and
d→∞
N
n1
e(1) → 0 and e(2) → 1 as d → ∞ if lim sup δ < − .
Nd→∞e(1) → 1 and e(2) → 0 as d → ∞ if lim inf δ >From Corollary 1, the DWD brings the strong inconsistency because of a huge bias
caused by heterogeneity of covariance matrices as well as sample imbalance. In order toovercome such diﬃculties, we proposed a bias-corrected DWD.

References[1] Marron, J.S., Todd, M.J., Ahn, J. (2007). Distance-weighted discrimination. Journalof the American Statistical Association, 102, 1267–1271.
[2] Qiao, X., Zhang, H. H., Liu, Y., Todd, M.J., Marron, J.S. (2010). Weighted distance weighted discrimination and its asymptotic properties. Journal of the American
Statistical Association, 105, 401–414.
[3] Vapnik, V. N. (2000). The Nature of Statistical Learning Theory (second ed.). New
York: Springer.
2順序カテゴリをもつ多元分割表における非同等性のモデリング
東京理科大学大学院 藤澤 健吾
東京理科大学 田畑 耕治

1. 正方分割表における一般化周辺非同等モデル
行変数 X と列変数 Y が順序のある同じ分類からなる r × r 正方分割表において, (i, j)
セル確率を πij (i = 1, . . . , r; j = 1, . . . , r) とする. 周辺同等 (MH) モデルは次のように表される (Stuart, 1955):
πi+ = π+i (i = 1, . . . , r),∑r∑r
ただし, πi+ = t=1 πit , π+i = s=1 πsi である. MH モデルは累積確率を用いて次のように表すこともできる:
FiX = FiY (i = 1, . . . , r − 1),∑i∑iただし, FiX = s=1 πs+ , FiY = s=1 π+s である. MH モデルがデータに当てはまらない
とき, MH モデルの拡張モデルを当てはめることに関心がある. 任意の k (k = 1, . . . , r −1)に対して, k 次一般化周辺非同等 (GMk ) モデルを次のように提案した:

G−1 (FiX ) =k−1∑

il ∆l + G−1 (FiY ) (i = 1, . . . , r − 1),

l=0

ここで G は G(−∞) = 0, G(∞) = 1 を満たす連続な単調増加関数である. 特に ∆0 =∆1 = · · · = ∆k−1 = 0 のときは MH モデルである. また, k = 1, G−1 (x) = log x/(1 − x)のときは周辺累積ロジットモデルである (McCullagh, 1977). さらに, k = 1, G−1 (x) =log(− log(1 − x)) のときは周辺累積補対数対数モデルである (Saigusa et al., 2018). その
他にも k = 1, G−1 (x) = Φ−1 (x) を考えることで周辺累積プロビットモデル等を構成することが可能である. ただし Φ(x) は標準正規分布の分布関数である.
任意の k (k = 1, . . . , r − 1) に対して, 周辺 k 次積率一致 (MEk ) モデルを次のように考
える:
r
r
∑
∑
l
i πi+ =
il π+i (l = 1, . . . , k).
i=1

i=1

このとき次の定理を得た:
定理 1 与えられた k (k = 1, . . . , r − 1) に対して, MH モデルが成立するための必要十分
条件は GMk モデルと MEk モデルの両方が成立することである.
この定理は Saigusa et al. (2018) 等の結果を含む.

2. 多元分割表における一般化周辺非同等モデル
Xt (t = 1, . . . , T ) が順序のある同じ分類からなる rT 分割表において, (i1 , . . . , iT ) セル
確率を P (X1 = i1 , . . . , XT = iT ) = πi1 ···iT (it = 1, . . . , r) とする. このとき rT 分割表における T 元周辺同等 (MH[T ]) モデルは次のように表される (Bhapkar and Darroch, 1990):
(1)

πi

(T )

= · · · = πi(i = 1, . . . , r),(t)

ただし πi = P (Xt = i) である. また累積確率を用いて次のように表すこともできる:

FiX1 = · · · = FiXT

(i = 1, . . . , r − 1),∑i(t)
ただし FiXt = s=1 πs である. MH[T ] モデルの当てはまりが悪い場合, MH[T ] モデル
より制約の弱いモデルに関心がある. 任意の k (k = 1, . . . , r − 1) に対して, T 元 k 次一般化周辺非同等 (GMk [T ]) モデルを次のように提案した:−1

G(FiXt ) =k−1∑

il ∆l + G−1 (FiX1 ) (i = 1, . . . , r − 1; t = 2, . . . , T ).
(t)

l=0

このモデルは Agresti (2002) 等のモデルを含む.
k (k = 1, . . . , r − 1) に対して, T 元周辺 k 次積率一致 (MEk [T ]) モデルを次のように考
える:
r
r
∑
∑
(T )
l (1)
(l = 1, . . . , k)il πi
i πi = · · · =
i=1

i=1

このとき次の定理が成立する:
定理 2 MH[T ] モデルが成立するための必要十分条件は GMk [T ] モデルと MEk [T ] モデルの両方が成立することである.
参考文献[1] Agresti, A. (2002). Categorical Data Analysis, 2nd edition. Wiley, New York.
[2] Bhapkar, V. P. and Darroch, J. N. (1990). Marginal symmetry and quasi symmetry ofgeneral order. Journal of Multivariate Analysis, 34, 173-184.
[3] McCullagh, P. (1977). A logistic model for paired comparisons with ordered categorical
data. Biometrika, 64, 449-453.
[4] Saigusa, Y., Maruyama, T., Tahata, K., and Tomizawa, S. (2018). Extended marginal
homogeneity model based on complementary log-log transform for square tables. International
Journal of Statistics and Probability, 7, 27-31.
[5] Stuart, A. (1955). A test for homogeneity of the marginal distributions in a two-way
classification. Biometrika, 42, 412-416.

High-dimensional classification by data transformation(データ変換を用いた高次元次判別分析について)
Aki Ishiia , Kazuyoshi Yatab , and Makoto Aoshimaba Department of Information Sciences, Tokyo University of Scienceb Institute of Mathematics, University of Tsukuba

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 that lim supp→∞ ||µi ||2 /p < ∞ for i = 1, 2, where || · || denotes the Euclidean
norm. Also, we assume that tr(Σi )/p ∈ (0, ∞) as p → ∞ for i = 1, 2. Here, for a function, f (·),“f (p) ∈ (0, ∞) as p → ∞” implies lim inf p→∞ f (p) > 0 and lim supp→∞ f (p) < ∞. We assumeni ≥ 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∑ni for i = 1, 2. Letz oj(i) = z j(i) − (z̄j(i) , ..., z̄j(i) )T , j = 1, ..., p; i = 1, 2, where z̄j(i) = n−1i
k=1 zjk(i) . We also
assume that()
P lim inf ||z o1(i) || ̸= 0 = 1 for i = 1, 2.
p→∞Also, we consider the following assumption as necessary:22(A-i) E(zqk(i)zsk(i)) = 1, E(zqk(i) zsk(i) ztk(i) ) = 0 and E(zqk(i) zsk(i) ztk(i) zuk(i) ) = 0 for all q ̸=s, t, u.

We note that (A-i) naturally holds when X i is Gaussian.
Let x0 be an observation vector of an individual∑
belonging to πi (i = 1, 2).∑We assume x0 and xij sii(xij − x̄ini )(xij −xij /ni and S ini = nj=1
are independent. We estimate µi and Σi by x̄ini = nj=1
Tx̄ini ) /(ni − 1). A typical classification rule is that one classifies an individual into π1 if{ det(S ) }2n2(x0 − x̄1n1 )T S −1
(x
−
x̄
)
−
log
< (x0 − x̄2n2 )T S −1
(1.1)
0
1n
1
1n1
2n2 (x0 − x̄2n2 ),det(S 1n1 )1and into π2 otherwise. However, the inverse matrix of S ini does not exist in the HDLSS context (p > ni ).
When Σ1 = Σ2 , Bickel and Levina (2004) considered the inverse matrix defined by only diagonal elements of the pooled sample covariance matrix. Yata and Aoshima (2012) considered using a ridge-typeinverse covariance matrix derived by the noise reduction (NR) methodology. When Σ1 ̸= Σ2 , Dudoit etal. (2002) considered using the inverse matrix defined by only diagonal elements of S ini . Aoshima andYata (2011,2014) considered substituting {tr(S ini )/p}I p for S ini in (1.1) by using the difference of ageometric representation of HDLSS data from each πi . Aoshima and Yata (2019) considered quadratic
classifiers in general and discussed asymptotic properties and optimality of the classifies under highdimensional settings. They showed that misclassification rates tend to zero as the dimension goes toinfinity. On the other hand, Chan and Hall (2009) and Aoshima and Yata (2014) considered distancebased classifiers and Aoshima and Yata (2014) gave the misclassification 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)tr(Σ2i )→ 0 as p → ∞ for i = 1, 2.

(1.2)

Recently, Aoshima and Yata (2018) considered the “strongly spiked eigenvalue (SSE) model” as follows:{ λ2
lim infp→∞1(i)tr(Σ2i )}> 0 for i = 1 or 2.

(1.3)

On the other hand, Aoshima and Yata (2018) called (1.2) the “non-strongly spiked eigenvalue (NSSE)
model”. For example, we consider a spiked model such as
λs(i) = as(i) pαs(i) (s = 1, ..., ti ) andλs(i) = cs(i) (s = ti + 1, ..., p)(1.4)

with positive and fixed constants, as(i) ’s, cs(i) ’s and αs(i) ’s, and a positive and fixed integer ti . Note that
(1.2) holds when α1(i) < 1/2 for i = 1, 2. On the other hand, (1.3) holds for the spiked model in (1.4)
with α1(i) ≥ 1/2.
As for the SSE model, Aoshima and Yata (2018) considered a classifier by using a data transformation from the SSE model to the NSSE model. They gave a consistency property of the classifier anddiscussed the asymptotic normality when p → ∞ and ni → ∞ (i = 1, 2).
In this talk, we created a new quadratic classifier under one of the SSE models. We emphasize thatone should construct a classification procedure by considering the eigenstructure of high-dimensional
data. We handled the divergence condition for p and ni ’s such as p → ∞ while ni ’s are fixed. We
introduced a quadratic classification procedure for high-dimensional data. We proposed a new quadratic
classifier under the SSE model. We showed that our classification procedure has a consistency property
for misclassification rates.

2

Robust Bayesian Modeling with
Synthetic Posterior∗
Shonosuke Sugasawa
Center for Spatial Information Science, The University of Tokyo

1

Introduction

The standard Bayesian method in which the posterior distribution is obtained as
the product of likelihood and prior would be able to work reasonably well when theassumed model is approximately correct. However, it is often the case that real datacontains outliers that the assumed model cannot explain. In this case, the standard
Bayesian method would break down in the sense that the posterior distribution does
not give reasonable Bayesian inference. In this project, we solve the diﬃculty ofthe standard Bayesian approaches by adopting synthetic posterior based on robust
divergence. In particular, we focus on Bayesian linear regression and dynamic linear
models and propose their robust alternatives.

2

Robust Bayesian Regression

Suppose we have independent observation (yi , xi ) for i = 1, . . . , n, where yi is acontinuous response and xi = (xi1 , . . . , xip )t is a p-dimensional vector of covariates.
We consider fitting a regression model yi = xti β+εi with ε ∼ N (0, σ 2 ). Let π(β, σ 2 ) be
prior distribution for the model parameters β and σ 2 . Then the standard posterior
distribution of the model parameters is given by the product of the prior and thelikelihood function, which is known to be sensitive to outliers, and it could produce
biased or ineﬃcient posterior inference.

∗

Joint work with Shintaro Hashimoto (Hiroshima University)1

To overcome the problem, we propose replacing the log-likelihood function withrobust alternatives. Specifically, we employ γ-divergence to define synthetic posteriorfor (β, σ 2 ). We develop an eﬃcient sampling strategy for the posterior distribution
via Bayesian bootstrap or nonparametric Bayesian updating, which generate posterior
samples by minimizing randomized loss function defined by the synthetic posterior.
The minimization problem can be eﬃciently carried out via MM-algorithm. In orderto address the case where the number of covariates is relatively large, we extended
the propose method to incorporated shrinkage priors on regression coeﬃcients.

3

Robust Dynamic Linear Modeling

Let (y1 , . . . , yt ) be a time series and we consider a simple dynamic linear model,
yt = θt + εt and θt = θt−1 + ut for t = 1, . . . , T , where εt ∼ N (0, σ 2 ), ut ∼ N (0, τ 2 ),and θt is a latent true signal. This model is useful for smoothing the observed timeseries yt to obtain latent trend (signal) θt . It is noted that θt is assumed to follow
random-walk process, in which the variance parameter τ 2 determines the smoothnessof θt . The unknown static parameters are σ 2 and τ 2 , which will be estimated byassigning prior distributions.
When εt contains outliers, the resulting signal estimator via the standard Bayesian
approach would be biased or ineﬃcient. To overcome the diﬃculty, we restrict the
information from observation yt by replacing the likelihood function ϕ(yt ; θt , σ 2 ) withrobust alternatives. Specifically, we here employ density power divergence to definesynthetic posterior distribution of latent signals. We developed a Gibbs samplingalgorithm for posterior computation.

2分割表における種々の対称モデルに関する尺度
東京理科大学大学院理工学研究科
横浜市立大学医学部
日本大学経済学部
東京理科大学理工学部

桃﨑 智隆 · 齋藤 健
三枝 祐輔
生亀 清貴
中川 智之 · 石井 晶 · 富澤 貞男

本講演では, 分割表における種々の対称モデルに関する尺度について 2 部構成で述べる.

第 1 部 正方分割表における部分対称と部分非対称からの隔たりを測る 2 次元尺度
正方 r × r 分割表を考え, (i, j) セル確率を pij (i = 1, · · · , r; j = 1, · · · , r) とする. {pij + pji > 0}を仮定する. Saigusa et al. (2016) は部分対称 (PS) モデルを提案した: 少なくとも 1 組の i ̸= j に∑∑対して, pij = pji . いま pcij = pij /(pij + pji ), p∗ij = pij /δ, δ =i̸=j pij とおく. Saigusa et al.
(2016) は PS モデルからの隔たりを測る尺度を提案した:

Φ(λ) =

r−1
∏

r
[
] ∗ ∗
∏
(λ) (pij +pji )ϕij(λ > −1),i=1 j=i+1

ただし,
(λ)

ϕij = 1 −

]
1[
(λ)
(λ)
c λ+1
c λ+1
.
H
,
H
=
1
−
(p
)
−
(p
)
ij
ji
ij
λ
2λ − 1 ij
λ2λ

本講演では, 新しいモデルと 2 つの尺度を提案する:(1) 部分完全非対称 (PAS) モデル:
pij = 0 (pji ̸= 0) または pji = 0 (pij ̸= 0) (少なくとも 1 組の i ̸= j)(2) PAS モデルからの隔たりを測る尺度:
τ (λ) =

r−1
∏

r
] ∗ ∗
[
∏
(λ) (pij +pji )1 − ϕij(λ > −1)i=1 j=i+1(3) 2 次元尺度:()⊤
Λ(λ) = Φ(λ) , τ (λ)(λ > −1)
Λ(λ) は次の性質を満たす:(i) Λ(λ) = (0, 0)⊤ ⇔ PS モデルと PAS モデルが同時に成立(ii) Λ(λ) = (0, 1)⊤ ⇔ 対称モデル (すべての i ̸= j に対して pij = pji ) が成立(iii) Λ(λ) = (1, 0)⊤ ⇔ 完全非対称モデル (すべての i ̸= j に対して pij = 0 (pji ̸= 0) または pji =0 (pij ̸= 0)) が成立
b (λ) は近似的に N (Λ(λ) , n−1 Σ) に従う. Λ(λ) に関するシミュレーション結果や, 信頼楕
推定尺度 Λ
円を用いた実データ解析については当日示す.

第 2 部 周辺点対称性からの隔たりを測る尺度
本節では，順序付き R × C 分割表における点対称性に着目する．ここで，R は行変数のカテゴリ
数，C は列変数のカテゴリ数を表し，各カテゴリは，行と列で順序付けがされているとする．第 1 部と同様に，(i, j) セル確率を pij (i = 1, · · · , R; i = 1, · · · , C) とし，i∗ = R+1−i, j ∗∗ = C +1−j と∑∑Rする．各 i, j に対し，行周辺確率を pi· = C
j=1 pij ，列周辺確率を p·j =i=1 pij とする．Tomizawa(1985) は，次の周辺点対称モデルを提案した：

pi· = pi∗ · (i = 1, . . . , R) かつ p·j = p·j ∗∗ (j = 1, . . . , C)

(1)

本講演では，周辺点対称モデル (1) からの隔たりを測る尺度を新たに提案する．{pi· + pi∗ · ̸= 0}と {p·j + p·j ∗∗ ̸= 0} を仮定する．

[ R2 ][ C2 ]∑∑
δ1 =(pi· + pi∗ · ), δ2 =(p·j + p·j ∗∗ )
i=1

j=1

とする．ここで，[x] は x を超えない最大の整数を表す．提案尺度は次のように定義される．

γMPS =
ここで，

δ1 γ1 + δ2 γ2
δ1 + δ2

[ R2 ]
(
)
∑
4
π
pi·
,
γ1 =
(qi· + qi∗ · ) θ1(i) −; θ1(i) = arccos  √
π
4
2
2
p +p
i=1



i∗ ·

i·


[ C2 ]
(
)
∑
∗∗
p·j
π
4

(q·j + q·j ∗∗ ) θ2(j) −; θ2(j) = arccos  √
γ2 =
π
4
2
p + p2 ∗∗
·j

i=1

であり，

p i∗ ·pi·, q i∗ · =
qi· =
δ1
δ1
p·j
p·j ∗∗q·j =, q·j ∗∗ =
δ2
δ2

·j

(
[ ])
R
i = 1, . . . ,
;
2
(
[ ])
C
j = 1, . . . ,2とする．提案尺度 γMPS は，以下の性質をもつ．

(i) −1 ≤ γMPS ≤ 1(ii) γMPS = 1 ⇔ γ1 = γ2 = 1 かつ γMPS = −1 ⇔ γ1 = γ2 = −1(iii) 周辺点対称モデル (1) が成立するならば，γMPS = 0 となる．
当日は，提案尺度の推定量における漸近正規性を述べ，近似信頼区間を構成し，実データへの適用例を紹介する．また，多元分割表に対する周辺点対称性について考え，新たな尺度を提案する．

参考文献
Saigusa, Y., Tahata, K., and Tomizawa, S. (2016). Measure of depature from partial symmetry
for square contingency tables. Journal of Mathematics and Statistics, 12:152–156.
Tomizawa, S. (1985). The decompositions for point symmetry models in two-way contingency
tables. Biometrical Journal, 27:895–905.

地上部と地下部のデータを統合させた推移行列モデル
島谷健一郎(統計数理研究所）、荒木希和子(立命館大生命科学部）
マルコフ推移行列モデルは、生物種の生活史のモデル化として、個体群生態学で積極的に活用されている。特に、行列要素に繁殖が加わることで、生物学固有の行列モデルが展開する。行列の行の和は 1 ではなく、固有値は生物個体群の増殖率を表すなど、推移行列モデルの数学と生物の間に、様々な対応が付けられる。
さらに、実データから行列要素を推定する統計モデルも数多く提唱されている。デー
タが多ければ単純な割り算しか必要ないが、野外の生物で大量のデータを得ることは困
難な場合が多い。、限られたデータから行列成分を推定する過程で統計モデルは必須となる。
植物は、大きく種子繁殖とクローン繁殖で拡大する。それぞれ行列モデルの中でモデル化されてきているが、実際の繁殖過程を大幅に簡略したモデルが主流である。とりわけクローン繁殖過程は、単なる年間平均クローン数程度に単純化されている。実際は、クローン繁殖する親株としない親株など、株の間に変異があるはずだが、モデル化以前に、どういう株がクローン繁殖するかといった基本的な野外データすら不十分な状況にある。
さらに、多年生の草本種には、地下茎でクローン繁殖する種も多い。その場合、どの地上部（株）がどの株の親株なのか、地上の観察ではわからない。遺伝子を調べること
でクローン識別はできるが、遺伝子情報だけでは親子関係の判別はできない。
地下を掘り起こし株と株のつながりを直接観察することで、親子関係は確かめられる。
ただし、1 度掘り起こせばその株は枯死する（または大きなダメージを被る）。1 度の破壊的調査で複雑なクローン繁殖の仕組みを定量的に知るには、掘り起こされた限られたデータから統計モデルによる推定が不可欠となる。
本研究では、クローナル繁殖動態について、(1) 数理モデルとしての定式化, (2) 実データによる検証法、を提唱することを主題に置いた。発表は、
1 地上部と地下部のフィールドデータ
2 クローナル繁殖動態のマルコフ推移行列モデルによる定式化
3 地上部と地下部のデータからの行列成分推定法
4. スズランのクローナル繁殖動態: クローン間の違いとその要因の順に行った。
2 では、以下のような構造の行列モデルを提唱した。
Cjh: age-position h, size class j の株のクローナル繁殖率
Stage: age-position (1-y-tip, 2-y-tip, non-tip, h) と size class (1 枚葉、2 枚葉 (<10cm,

10-15cm, 15-20cm, 20cm<, i, j)

Cj1-y-tip推移
非先端

Cj2-y-tip

1 ‒ Cj2-y-tip

Cnon-tip1 ‒ Cj1-y-tip

2年目以降
先端

(分枝)

Cj2-y-tip

1年目先端

Cj1-y-tip

Cnon-tip (分枝)繁殖さらに、以下の記法を加えて、行列成分は以下のようになる。
Ni; サイズクラス i の新規クローン割合
Sj: サイズクラス j の生残率
Tij:サイズクラス j が生残したときクラス j に推移する確率
E: 平均娘数（高々1 シーズンに 3 クローン）Year t (当年）size class i

繁殖によるage-position変化 (C ), 生残率 (S ) , 推移確率（T )Year t - 1 (前年), size class j
非先端
1年目先端 1-y-tip
2年目先端 2-y-tip
non-tip

繁殖による増殖確率 (C ), 新規サイズ分布 (N j ), 新規娘数 (E )Year t - 1 (前年), size class j
1年目先端 1-y-tip

2年目先端 2-y-tip

非先端 non-tip

EN i C j 1-y-tip

EN i C j 2-y-tip

EN i C j non-tip

1年目先端
1-y-tip

0

0

0

1年目先端
1-y-tip

2年目先端
2-y-tip

T ij S j (1 − C j 1-y-tip )T ij S j (1 − C j 2-y-tip )0+ 2年目先端2-y-tip000

非先端 nontip

0

0

0

非先端 nontip

T ij S j C j

1-y-tip

T ij S j C j

2-y-tip

T ij S j

実データから各成分をベイズ推定し、得られた行列の固有値の分布を求めたところ、
以下のように 3 つのクローン間で異なる分布を示す結果となった。
1400

1010固有値
1020固有値

1200

1031固有値

1000
800
600
400
200
0
1.2-1.205
1.195-1.2
1.19-1.195
1.185-1.19
1.18-1.185
1.175-1.18
1.17-1.175
1.165-1.17
1.16-1.165
1.155-1.16
1.15-1.155
1.145-1.15
1.14-1.145
1.135-1.14
1.13-1.135
1.125-1.13
1.12-1.125
1.115-1.12
1.11-1.115
1.105-1.11
1.1-1.105
1.095-1.1
1.09-1.095
1.085-1.09
1.08-1.085
1.075-1.08
1.07-1.075
1.065-1.07
1.06-1.065
1.055-1.06
1.05-1.055
1.045-1.05
1.04-1.045
1.035-1.04
1.03-1.035
1.025-1.03
1.02-1.025
1.015-1.02
1.01-1.015
1.005-1.01
1-1.005
0.995-1
0.99-0.995
0.985-0.99
0.98-0.985
0.975-0.98
0.97-0.975
0.965-0.97
0.96-0.965
0.955-0.96
0.95-0.955
0.945-0.95
0.94-0.945
0.935-0.94
0.93-0.935
0.925-0.93
0.92-0.925

また、いくつかのモデルを比較・対照させることにより、クローン間の違いは、生残・推移（成長）よりクローナル繁殖率に起因すると思われた。