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

大規模複雑データの理論と方法論：最前線の動向と新たな展開
科学研究費補助金 基盤研究 (A) 20H00576「大規模複雑データの理論と方法論の革新的展開」（研究代表者: 青嶋誠），学術研究助成基金助成金 挑戦的研究 (萌芽) 19K22837「超高次元データによる個別化モデリン
グへの挑戦」（研究代表者: 青嶋誠）によるシンポジウムを下記のように催しますので，ご案内申し上げます．

誠(筑波大学)

矢田和善(筑波大学)

石井(東京理科大学)

青嶋晶

記
日時：2020 年 12 月 4 日 (金)〜5 日 (土)場所：つくば国際会議場 大会議室 102 ( https://www.epochal.or.jp/ja/ )〒305-0032 茨城県つくば市竹園 2-20-3

プログラム
12 月 4 日 (金)10:00〜10:10

開会

10:10〜10:50

中山

優吾（京都大・情報学研究科）

矢田

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

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

高次元カーネル主成分分析に基づく異常値の検出

11:00〜11:40

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

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

順序カテゴリをもつ二元分割表における独立性の必要十分条件

11:50〜12:30

庄建倉（統計数理研究所）

Model Improvement Through Residual Analysis
12:30〜14:00

昼食

14:00〜14:40

永井勇（中京大・教養教育研究院）

rank 落ちを回避する新たな推定法とその応用
14:50〜15:30

柳本武美（統計数理研究所）

ベイズ型対数尤度に基づくモデルの信用集合

15:45〜16:25

竹内努（名古屋大・理学研究科）

Analysis of Spatially Resolved Galaxy Spectra as a High-Dimensional Low-Sample Size
Data Problem
16:35〜17:15

竹内努（名古屋大・理学研究科）

Topological Data Analysis of Galaxy Distribution: Detection of the Birth Cry of
Cosmic Structures

12 月 5 日 (土)9:00〜9:40

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

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

Sparse PCA by the noise-reduction methodology
9:50〜10:30 沈 迅（東京農工大・工学部）

庄 建倉（統計数理研究所）

Decimating Nonlinear Response in State Space Models for Battery Capacity Estimation
10:40〜11:20

地嵜頌子（大阪工業大・情報科学部）

宮本暢子

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

藤原良叔（筑波大学）

二部グラフ構造を持つ処理集合の実験計画とその深層学習への応用

11:35〜12:15

大久保祐作（統計数理研究所）

マクロ進化の大規模種間比較研究における枝特異的方向性淘汰モデルの提案

12:25〜13:05

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

２次林再生における林冠個体の 15 年間変遷を観る統計モデル

13:05〜13:10

閉会高次元カーネル主成分分析に基づく異常値の検出

京都大学・情報学研究科
筑波大学・数理物質系
筑波大学・数理物質系

1

中山 優吾（Yugo Nakayama）
矢田 和善（Kazuyoshi Yata）
青嶋 誠（Makoto Aoshima）

はじめに

本講演では，高次元小標本データに対する異常値の検出を考えた．高次元データの異常値の検出について，マハラノビス距離に基づく異常値の検出が盛んに研究されている．しかしながら，高次元小標本において標本共分散行列の逆行列が存在しないため，マハラノビス距離を単純に用いることができない．本講演では，高次元カーネル PCA を用いたクラスタリングを応用して，高次元小標本データに対して異常値検出を考えた．

2

主成分スコアにおける高次元一致性

2 クラスの高次元クラスタリングを考えることで異常値の検出を提案した．2 つの d 次元分布を Π1 , Π2と名付け，それぞれ平均 µ1 , µ2 と，共分散行列 Σ1 , Σ2 をもつと仮定する．ここで，tr(Σ1 ) ≤ tr(Σ2 )とする．それらの混合分布から n (≥ 2) 個のデータを無作為に抽出し，データ行列を X = [x1 , ..., xn ]とする．ただし，ni = #{j|xj ∈ Πi , j = 1, ..., n} とおく．ここで，#A は集合 A の要素の個数とし，ni ≥ 1 とする．
本節では，線形カーネル
k(xi , xj ) = xTi xj
を用いたカーネル PCA を考える．I n を n 次の単位行列，1n = (1, ..., 1)T として，P n = I n −n−1 1n 1Tnとする．カーネル関数 k(xj , xj ′ ) を (j, j ′ ) 成分にもつグラム行列 K を考え，その中心化グラム行列を
K 0 = P n KP n とし，K 0 の固有値分解を

K0 =

n−1
∑

λ̂i ûi ûTi

(ûi = (ûi1 , ..., ûin )T ,∥ûi ∥2 = 1)i=1とする．xj の (基準化した) 第 i 主成分スコアを sij =√nûij とおく．次を仮定する．

(A-i) Var(∥x − µi ∥2 |x ∈ Πi ) = O{tr(Σ2i )} as d → ∞ for i = 1, 2.
Πi が正規分布ならば，Var(∥x − µi ∥2 |x ∈ Πi ) = 2tr(Σ2i ) となるので，(A-i) は正規性の条件を緩めた仮定であることに注意する．さらに，∆µ = ∥µ1 − µ2 ∥2 ，∆Σ = |tr(Σ1 ) − tr(Σ2 )| とおく．次を仮定する．
(A-ii) tr(Σ2i )/∆2µ = o(1) as d → ∞ for i = 1, 2.
そのとき，次の定理が成り立つ．

1

定理 1 (Nakayama et al. (2020)). (A-i) と (A-ii) を仮定する．さらに，次を仮定する．

(A-iii) lim supd→∞∆Σ< 1 when n2 ≥ 2.
n1 ∆µd → ∞ のとき次が成り立つ．

{√
n /n + oP (1)
√2 1
s1j =
− n1 /n2 + oP (1)when j = 1, ..., n1 ,when j = n1 + 1, ..., n.

(1)

注意 1. Yata and Aoshima (2019) は 定理 1 と同様の結果を ∆Σ /∆µ = o(1) のもとで与えた．(A-iii)は ∆Σ /∆µ = o(1) よりも緩い条件であることに注意する．n2 = 1 のとき，(A-i) と (A-ii) のもとで
(1) が成り立つ．
注意 2. もし tr(Σ1 ) > tr(Σ2 ) なら，(A-i), (A-ii) と次の条件のもとで (1) が成り立つ

lim sup ∆Σ /(n2 ∆µ ) < 1 when n1 ≥ 2.
d→∞

定理 1 から (A-iii) のもとで，第 1 主成分スコアの符号で高次元データを分類できることがわかる．
しかし ∆µ /∆Σ が小さい場合，線形カーネルでは分類することはできず，非線形カーネルを使う必要がある．例えば，ガウシアンカーネル

k(xi , xj ) = exp(−∥xi − xj ∥2 /γ) (γ > 0)を用いたカーネル PCA では，µ1 = µ2 (∆µ = 0) のときでさえも共分散行列間の距離によって一致性(1) を与えることができる．Nakayama et al. (2020) は高次元におけるカーネル主成分スコアの漸近的性質を導出し，平均ベクトル間の距離 ∆µ だけでなく，共分散行列間の距離 ∆Σ にもよって (1) のような主成分スコアの高次元一致性をもつことを示した．
本講演では，主成分スコアの一致性 (1) を異常値検出に応用した．例えば，異常値が 1 つと想定できる場合は，n1 = 1 として考えることで (1) より異常値の検出が可能である．当日は，一般のカーネ
ル関数の場合と異常値が複数ある場合についても言及し，主成分スコアの符号に基づき異常値検出が可能であることを理論的に示し，その性能を数値実験と実データ解析を用いて確認した．

References[1] Nakayama, Y., Yata, K. and Aoshima, M. (2020). Clustering by principal component analysis
with Gaussian kernel in high-dimension, low-sample-size settings, submitted.
[2] Yata, K. and Aoshima, M. (2019). Geometric consistency of principal component scores forhigh-dimensional mixture models and its application. Scand. J. Stat., 47: 899–921.

2

順序カテゴリをもつ二元分割表における独立性の必要十分条件
東京理科大学大学院 藤澤 健吾
東京理科大学 田畑 耕治

1. 二元分割表におけるモデル
順序カテゴリをもつ r × c 分割表において, (i, j) セル確率を πij (i = 1, . . . , r; j = 1, . . . , c)とする. 独立 (I) モデルは次のように表される:

πij = πi+ π+j (i = 1, . . . , r; j = 1, . . . , c),∑c∑r
ただし, πi+ = t=1 πit , π+j = s=1 πsj である.
I モデルの当てはまりが悪い場合, I モデルよりも制約の弱いモデルを当てはめることに関心がある. Goodman (1979) は一様連関 (U) モデルを次のように提案した:
πij = µαi βj θij(i = 1, . . . , r; j = 1, . . . , c).

U モデルは θ = 1 のとき, I モデルである. また Agresti (1983) は線形-線形連関 (LL) モ
デルを次のように提案した:
πij = µαi βj θui vj(i = 1, . . . , r; j = 1, . . . , c),

ただし {ui } は u1 < · · · < ur を満たす行カテゴリに割り振られた既知のスコアであり,{vj } は v1 < · · · < vc を満たす列カテゴリに割り振られた既知のスコアである.
p = (pij ), q = (qij ) を離散確率分布とする. 確率分布 p と q に対する f -ダイバージェンスは次のように表される (Csiszár and Shields, 2004):
( )
∑∑
pij
I(p||q) =
qij f
,
qij
i
j
ただし f は f (1) = 0, f (0) = limt→0 f (t), 0 · f (0/0) = 0, 0 · f (a/0) = a limt→∞ (f (t)/t) を満たす (0, ∞) 上の凸関数である.
f を 2 回微分可能な狭義凸関数とする. Kateri and Papaioannou (1994) は f -ダイバー
ジェンスに基づく線形-線形連関 (LL[f ]) モデルを次のように提案した:

πij = πi+ π+j F −1 (αi + βj + θui vj ) (i = 1, . . . , r; j = 1, . . . , c),
ただし, F (t) = f ′ (t) である. f (t) = t log t のとき, LL[f ] モデルは LL モデルである.

3. 独立性の必要十分条件
行スコアに対する確率変数を U , 列スコアに対する確率変数を V とする. このとき,∑∑E(U ) =
i ui πi+ , E(V ) =j vj π+j である. Tomizawa et al. (2008) は独立性の必要十分条件に関する次の定理を与えた:

定理 1 独立性が成立するための必要十分条件は LL モデルと Cov(U, V ) = 0 が同時に成立することである.
本報告では LL モデルの一般化である LL[f ] モデルを用いて次の定理を与えた:
定理 2 独立性が成立するための必要十分条件は LL[f ] モデルと Cov(U, V ) = 0 が同時に成立することである.

3. モデルの直交性
nij を (i, j) セル観測度数, mij を (i, j) セル期待度数とする. また, {nij } は総観測度数∑ ∑n = i j nij の多項分布に従うことを仮定する. モデル M に対する尤度比カイ二乗検
定統計量は次のように表される:
(
)
c
r ∑
∑
nij
2
G (M ) = 2
nij log,m̂ij
i=1 j=1
ただし, m̂ij は mij の最尤推定量である.
行スコア {ui } と列スコア {vj } にE(U ) = E(V ) = 0を仮定する. このとき次の定理を与えた:
定理 3 独立モデルのもとで, G2 (I) は G2 (LL[f ]) と G2 (Cov(U, V ) = 0) の和に漸近的に同等である.
参考文献[1] Agresti, A. (1983). A survey of strategies for modeling cross-classiﬁcations having
ordinal variables. Journal of the American Statistical Association, 78, 184–198.
[2] Csiszár, I. and Shields, P. C. (2004). Information Theory and Statistics: A tutorial.
Now Publishers, Hanover.
[3] Goodman, L. A. (1979). Simple models for the analysis of association in crossclassiﬁcations having ordered categories. Journal of the American Statistical Association,
74, 537–552.
[4] Kateri, M. and Papaioannou, T. (1994). f -divergence association models. International Journal of Mathematical and Statistical Sciences, 3, 179–203.
[5] Tomizawa, S., Miyamoto, N., and Sakurai, M. (2008). Decomposition of independence
model and separability of its test statistic for two-way contingency tables with orderedcategories. Advances and Applications in Statistics, 8, 209–218.

Model improvement through residual analysis
Jiancang Zhuang
The Institute of Statistical Mathematics, Tokyo, Japan

Statistical modelling has been increasingly and widely used in many research fields. Associated statisticalinference techniques, such as model selection, specification, and evaluation, parameter estimation, and testing
goodness-of-fit have been developed. Among them residual analysis can be used to evaluate a model andassess the goodness-of-fit without proposing unnecessary new models whose implementation may include
heavy programming and computation tasks (Zhuang, 2006; Baddeley et al., 2005). Residual analysis is also
helpful for understanding the advantages of a model even if its overall performance is not ideal.
More recently, general principles of first-order residuals for spatial point processes (Papangelou intensitybased models), temporal/spatiotemporal point processes and moment intensity-based models were proposedby Baddeley et al. (2005), Zhuang (2006), and Guan et al. (2008), respectively. Zhuang (2015) summarized
the general principles for constructing innovation based residuals. For hidden Markov models, Buckby et al.
(2020) developed a similar innovation based residuals by extending the framework of residual analysis in
Zhuang (2015).

1

Residual analysis for point processes

A point process N is a random counting measure on Rd satisfying that N (A) for every bounded set A ⊂ Rd .
In this article, all point processes are assumed to be simple; that is, each point process can be represented bya sum of unit Dirac measures at distinct locations.
For a point process N on E, its first-order moment intensity, provided that it exists, is defined on E byµ(x) = limδ↓0

N (B(x, δ)).
`(B(x, δ))(1)where ` is the Lebesgue measure, and B(x, δ) represents the ball centered at x and with radius δ. The following
fact holds.
Conditional intensities When a temporal axis is embedded in Rd , i.e., the other d − 1 dimension are
regarded as spatial directions or marks. The conditional intensity function λ(t, x), t ∈ R and x ∈ Rd−1 , isdefined as
E [N ([t, t + δt ) × B(x, δx )) | Ht ]λ(t, x) = lim(2)
δt ↓0,δx ↓0
δt `(B(x, δx ))where Ht represent the observational history of N up to t but not including t.
Papangelou intensity The Papangelou intensity is defined by, if it exists,
λp (x) = lim

δx ↓0

E [N (B(x, δx )) | E(B(x, δx ))]δt `(B(x, δx ))(3)

where E(B(x, δx )) represent the observation on N outside of the close ball B(x, δx ).
The first-order residuals for moment intensity models, conditional intensity models, and Papangelou intensity models are defined by
Z hi
Rm (B, f ) =fˆ(x) N (dx) − fˆ(x) µ̂(x) `(dx) ,
(4)
B

Z h
i
Rh (B, f ) =fˆ(x) N (dx) − fˆ(x) λ̂(x) `(dx) ,Z Bhi
Rp (B, f ) =fˆ(x) N (dx) − fˆ(x) λ̂p (x) `(dx) ,

(5)
(6)

B

respectively, where B is a Borel set, and f is a deterministic function in (4), a predictable function in (5), andan exvisible function in (6). If f also depends on the estimated parameters, then fˆ is obtained by substituting
the estimated parameters in f . In practice, it is reasonable to assume that
R{m,h,b} (B, f ) ≈ 0,
when the estimated model is a good approximation of the true model.
1

(7)

2

Residual analysis for Hidden Markov models

Give an HMM whose the observations are realised from the m members of a family of distributions, f (xt |St = i), with i = 1, · · · , m, where xt is the observed value at time t = 1, · · · , T , St is the state of the Markov
chain at time t and m is the finite number of states with a stationary transitive matrix
Γij = P (St = j | St−1 = i).
In an HMM, the forecasted distribution of Xt , given the previous observations is,m PmX(k) Γ f (Xt | St = i)k=1 αt−1
Pm ki,
f (xt | Ht ) = f (xt | X1 , · · · , Xt−1 ) =k=1 αt−1 (k)
i=1

(8)

(9)

where αt (k) is the forward probability defined by
αt (k) = P (X1 , · · · , Xt , St = k)

(10)

Using the expectation calculated above, the predictive residual with a predictable function h(t, x), where
h is a predictable function of t for every fixed x, becomes,
Z
n
X
Rn =
(h(t, Xt ) − h(t, xt )f (xt | Ht ) dxt .
(11)
t=1

For the exvisive residuals, we consider the forecasted distribution of Xt given all of the other observations
rather than observations prior to time t,m PmXt−1 (k) Γki f (xt | St = i) βt (i)k=1 αPf (xt | Et ) =(12)m
k=1 αt−1 (k) Γki βt (i)i=1where βt (i) = P (Xt+1 , · · · , XT | St = i) is the backward probability. Given a function h(t, x) which isexvisible as a function of t for every x, the exvisible residual becomes

Z
n 
X
e
Rn =
h(t, Xt ) − h(t, xt ) f (xt | Et ) dxt .
(13)
t=1

It is easy to see E[Rn ] = 0 and E[Rnp ] = 0.

3

Summary

Based on the theory of residual analysis for point processes and hidden Markov models, the stochastic reconstruction technique is illustrated with two examples: stochastic reconstruction of a Hawkes point processmodel for the robbery related violence in Castellon, Spain, and residual analysis to the results from fitting azero-inflation HMM to the tremor data in the Kii region, Japan. Residual analysis helps us to ?nd a proper
direction in which to improve on the model formulation, in a data-driven manner.

References
Baddeley, A., Turner, R., Møller, J., and Hazelton, M. (2005). Residual analysis for spatial point processes(with discussion). Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67(5):617–
666.
Buckby, J., Wang, T., Zhuang, J., and Obara, K. (2020). Model checking for hidden Markov models. Journalof Computational and Graphical Statistics, 0(0):1–16.
Guan, Y., Waagepetersen, R., and Beale, C. M. (2008). Second-order analysis of inhomogeneous spatial
point processes with proportional intensity functions. Journal of the American Statistical Association,
103(482):769–777.
Zhuang, J. (2006). Second-order residual analysis of spatiotemporal point processes and applications in model
evaluation. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(4):635–653.
Zhuang, J. (2015). Weighted likelihood estimators for point processes. Spatial Statistics, 14(B):166–178.
2

ベイズ型対数尤度に基づくモデルの信用集合
柳 本 武 美 統計数理研究所

 序

 

 

次元多項分布 Ô は自然母数を と置いた曲指数型分布と見なせる。多項分布モデルに還元させるためには、自然母数と結びつくようにする。そのためにモデルの尤度ではなくモデルの対数尤度を定義する必要がある。異なる
個の事前分布をもつベイズモデル対して、そのベイズ型対数尤度を対数尤度の事後平均により定義する。 個のベイズモデルは、分類問
題の個のクラスと定式化してそのクラスに事前確率を仮定する。この定式化により、事後確率が誘導されモデルの信用集合が得られる。信用集合の導入により得られたモデルの評価が可能になる。例えば、仮説検定とモデル選択の手続きが統一的に得られる。単にモデルの最適な分類を超えた幅広い適用が期待される。

 ベイズ型対数尤度
標本密度として指数型分布を仮定する。


Ü





    

 Ü
  

定して、それらを



 







  

Ü



 

  以降では、標本密度は固定して議論する。異なる  個の事前密度を仮 とする。これから Ü       

 







 



個のベイズモデル

 



 

が得られる。ベイズ型対数尤度は、対数尤度の事後平均 ÔÓ×Ø  Ü  が基本となる。しかし、このままでは標本 Ü の二重使用になっているのでその偏りを補正する必要がある。
その補正のためには ÔÓ×Ø  Ý   を評価する。簡単のために漸近近似が使われるが、有限標本での妥当性に問題が生じる。。そこでこの事後平均が最適な予測子と関係することに注意する。予測子を Ý  ÜÔÓ×Ø  Ý   で定義すると、この予測子は  Ý ÜÝ   の事後平
均を最小にする意味で最適な予測子である。更にその予測子は の事後平均 を  して最適な予測子が得られる。つまり  Ý  ÜÝ   となる  !"! ##$ %&'(。これからモデルのベイズ型対数尤度  
  を

 



 











 

ÔÓ×Ø 

 




 Ü


    Ý Ü  Ý 














 

)

で定義する。この尤度は次の別表現が可能である。
 Ü Ü Ý Ü     ÔÓ×Ø 







 



 モデルの事前分布と信用集合
尤度ではなくて対数尤度を扱った理由のは次元多項分布 Ô の構造にある。多項分布は自然母数を  Ô  ½
   と置いた正値測度空間の部分空間と見なせる。解析的には曲指数分布族と同じように扱える。そこで、ベイズ型対数尤度を用いて  










 



と置いて、多項分布の推定問題に還元させる。母数 Ô は  次元標準単錐  ½ の要素である。その単錐の頂点全体を  ½ と書くと、 個の各々の頂点は分類問題におけるクラスに対応する。こ ½ とする。事前確率の問題ではモデルに対応する。そこでモデルの事前確率を   Õ 
が * + 選択になるが、問題に応じて自由に選ぶことが出来る。これから事後確率
は









Ü

 
















 
 






が得られる。
事後確率から信用水準、信用区間が導かれる。本節のような手続きは無用に感じるかも知れない。
ごく常識的は




,-





 






として  を求めれば十分と考えるからである。しかし実際には   Ü の値の大きさが分類結果のに近いと分類の意味は弱い。また、信用水準
解釈に大きな影響を及ぼす。例えば、この値がを与えておいて、適当な に対して .' 信用集合は    Ü

が







¾



Ü







 





を満たすように選ぶことで定義される。実用上の分類問題で必要になる情報を提供する。。

 適用
  モデル選択
モデル選択は、文字通りのモデルの選択と選択したモデルの下での母数の推定と推定した母数を した予測子の構成の二面がある。しかし、その理論では両者は一体的である。伝統的には /(0と表され多くの手法、 ,(0 1(0 (0 2,(0 などが提案されている。大雑把近似では  3/(0がベイズ型対数尤度に相当する。尤度法におけるモデル選択では  ではなく母数の集合  により定義される。  は


 Ü




 



   




  



Ü   

 

 

と表される。モデルは 4 Ü    である。事後確率に相当する量の研究が行われ 5-!
,* -" ## はそのレビューを与えている。モデルの事前分布は、モデルの重みとされる。定義の理論体系が曖昧である。
提案法の特徴は 関連した予測子が標本サイズ有限の下で最適性がある、 ,(0 1(0 (0 と比較して 漸近的近似は援用しない) 観察予測子Ü Ü ではなく対数尤度を基本とする
6 モデルを与えると一意に定まる、 (0 と比較して 仮説検定
節でも議論したようにベイズ理論における信頼区間に対応する信用区間を定義した。従って、モデルの事前確率を選ぶことにより、仮説検定を説明できる、例えば - # 7 & *。その結果、モデル選択の理論で強調されている仮説検定との乖離がなくなる。既存の研究では、ベイズ因子とか /(0 によるモデルの選択が仮説検定よりも優れていることを主張されていた。
本研究では、優劣を論じないで両者を統一的に俯瞰する。
)

Analysis of Spatially Resolved Galaxy Spectra as aHigh-Dimensional Low-Sample Size Data Problem(高次元小標本データ問題としての銀河の分光マップの解析)
竹内 努 1,2, 河野 海 1, 中西 康一郎 3, 矢田 和善, 青嶋 誠 4, 石井 晶 5, 江頭 健斗 6
1. 名古屋大学理学研究科素粒子宇宙物理学専攻,
2. 統計数理研究所統計的機械学習研究センター
3. 国立天文台 ALMA プロジェクト
4. 筑波大学数理物質系
5. 東京理科大学理工学部情報科学科
6. 筑波大学 理工情報生命学術院
銀河とは, 星と星間物質(ガスと星間塵), 暗黒物質からなる巨大な天体である. 銀河は現在観測可能な宇宙に数千億個存在しており, 我々の目に見える波長(可視光線)での宇宙の姿を形作っている. しかし宇宙誕生当時, 物質分布はほぼ一様であった. 銀河は平均からわずかに密度が高い領域が重力によって成長し, 合体成長を経て現在の姿へと進化してきたのである. 銀河進
化は周囲の銀河の密度やガス密度など, 銀河の置かれた環境にも強く依存する極めて複雑な過程である. そして銀河の形成と進化は 138 億年の宇宙進化の歴史の中でも非常に重要な現象のひとつであると考えられている.
銀河の進化を駆動する重要な内的物理過程が星形成である. 星形成は星間物質が重力によって凝集し, 電離ガスから原子ガス, そして分子ガスへと冷却することで分子雲, さらに高密度のガス塊を形成し, 核融合が始まるという過程である. 星になる母体が分子雲であることから, 分子雲の物理状態を解明することが銀河進化の重要な課題となっている. 観
測的にも星形成率がガス密度の冪乗に比例するとい う 経 験 則 (Kennicutt-Schmidt 則 : Schmidt 1959;

図 1: 電離ガス(プラズマ)から中性ガス雲、分子雲への凝集。

Kennicutt 1989, 1998)が知られており, 星間分子の素過程の解明はこの経験則を物理的に説明するための手がかりとしても重視されている.
実験室実験の不可能な遠方宇宙を対象とする天文学において, 天体の分子の物理情報を得る唯一の手段が分光観測である. 分子の特徴的量子状態に対応する電磁波放射(輝線放射)は
主として電波に現れるため, 分子雲の物理は電波天文学で発展してきた. 分子の電波輝線を観
測 で き る , 現 在 最 も 大 規 模 か つ 高 性 能 の 電 波 観 測 装 置 が Atacama Large
Millimeter/Submillimeter Array (ALMA)である. ALMA を始めとする現在の天体分光観測装置は,天体の星間物質が含む原子, 分子, イオン, そしてダスト(星間塵)について圧倒的な質・量の情報を提供する. しかし詳細な分光観測は一般的に非常に時間がかかり, 天体をマッピングして独立な観測点を多数得ることは容易ではない. その結果, データの持つ波長方向の次元 d に対してサンプル数 n が非常に小さい(n ≪ d)という状況が生じる. 伝統的な天文学では, このような状態は不良設定問題とされ, その解決は諦めざるを得ないとされてきた. 物理的考察のためにはデータ d の情報を大幅に捨て, n より小さくして解析する他ないと考えられていたのである.
当然このような無駄は避け, 情報量をフルに活かせる方法が望ましい. n ≪ d となるデータは高次元小標本(HDLSS)データと呼ばれる. 天文学以外の分野, たとえばゲノム解析では遺伝子の標本数 n ~ 100 に対し塩基配列の次元数が d ~ 105 といった問題は特殊ではない. HDLSS の問題を解決するためにここ 10 年で発展した統計学の方法が高次元統計学であり, 現在も続々と新しい知見が得られている(e.g., Aoshima 2018). 本研究では, 天体の各点を ALMA によって分光観測した探査(マッピング)型のデータに対して高次元統計学の方法を適用した新しい解析を行っている. 上記のように, 分光マッピングには大変な観測的コストがかかり, ALMA の場合現状では最大級でも n ~ 200 程度の観測点が得られる程度であるのに対して波長(振動数)分散方向に
は d ~ 2000 の情報がある. よって, ALMA の分光マッピングは典型的な HDLSS データである. 従来の天文学では, 分光データから物理の第一原理から有用であるとわかっているいくつかの輝線のみを抽出し, それを用いて分類を行っていた. しかし, ALMA の観測によって, たとえ隣り合った星間分子雲の領域であっても分子輝線は大幅に異なった様相を示すことが判明し(図 2), 古典
的方法の限界が明らかになっている. 巨大すぎる情報量のため, 分子雲の進化と星形成の関係が全く見えてきていないのである. 高次元統計学の方法によって, かつて絶望的と考えられていた, 星間物質の輝線の全ての情報を活かした分類が可能になると期待される. もちろん, 天文学
での分光データに対する高次元統計学の適用は前例がなく, 天体分光固有のデータの問題
あり, その可能性は未知数の部分も大きい. 本講演では問題意識と目標, そして研究の現状について述べる.

図 2: ALMA による爆発的星形成銀河 NGC253 中心部の星形成領域の分光マッピングの例(Ando et al. 2017より改変). 右パネルが観測領域のマップ、左パネルが観測点に対応する分子輝線スペクトルデータを示す。

Topological Data Analysis of Galaxy Distribution:
Detection of the Birth Cry of Cosmic Structures(銀河分布の位相的データ解析: 宇宙構造の産声の検出)
竹内 努 1,2, 河野 海 1, Suchetha COORAY1, 西澤 淳 3,1, 村上 広椰 1
1. 名古屋大学理学研究科素粒子宇宙物理学専攻,
2. 統計数理研究所統計的機械学習研究センター
3. 名古屋大学高等研究院
銀河とは, 星と星間物質(ガスと星間塵), 暗黒物質からなる巨大な天体である. 銀河は現在観測可能な宇宙に数千億個存在しており, 我々の目に見える波長(可視光線)での宇宙の姿を形作っている. しかし宇宙誕生当時, 物質分布はほぼ一様であった. 銀河は平均からわずかに密度が高い領域が重力によって成長し, 合体成長を経て現在の姿へと進化してきたのである(図 1). 銀河進化は周囲の銀河の密度やガス密度など, 銀河の置かれた環境にも強く依存する極めて複雑な過程である. そして銀河の 図 1: 一様な物質分布から銀河が形成さ形成と進化は 138 億年の宇宙進化の歴史の中でも非 れ、進化する様子.
常に重要な現象のひとつであると考えられている.
我々は 130 数億年にわたる銀河の形成・進化という複雑な物理現象を, これまでの宇宙物 理 学の 方法 とはまったく異なる, データ科学の最新手法である位相的データ解析(topological data analysis: TDA)に基づくアプローチによって新たな角度から解明することを試みている. 本研究は特に, 銀河形成以前の通常物質(宇宙物理学ではバリオンと呼ぶ)が
光 と の 相 互 作 用 に よ っ て 生 じ る 粗 密 波 で あ る 「 バ リ オ ン 音 響 振 動 」 (baryon acoustic
oscillation: BAO)が銀河分布に残す痕跡について扱う. バリオンと光の相互作用が切れる宇
宙の晴れ上がり(脱結合)時に物質が放射した熱放射は宇宙膨張によって波長が引き伸ばされ, 現在は電波(ミリ波)の波長で宇宙マイクロ波
背景放射(cosmic microwave background: CMB)として観測される. BAO は CMB のゆらぎに特徴的なスケールとして現れ, ゆらぎをフーリエ変換
したパワースペクトルにおいて最も顕著なピー
クとして観測される(図 2). このスケールは共動距離で観測した際に, どの赤方偏移においてもおよそ 148 Mpc (1 Mpc = 3.08×1024 cm)の一定値
図 1: Planck 衛星が観測した CMB のパワースをとる. この性質から, BAO シグナルの観測的ペクトル. 最も顕著なピークが BAO である.

スケールは宇宙論パラメータの推定の最も信頼できる指標のひとつとなっている.
従来, BAO の検出には銀河分布の 2 点相関関数やパワースペクトルが用いられてきたが,BAO シグナルは銀河の重力的クラスタリングによって生じるシグナルに対して 1/100 程
度と微弱であるため, 統計的に有意な結果を得るためには大規模かつサンプリングレートの高い銀河探査が必要であり, 現存する最大の銀河探査であるスローンデジタルスカイサーベイ(Sloan Digital Sky Survey: SDSS)によって初めてそのシグナルが検出された. この規模の探査は決して容易ではなく, 独立な, そしてより有効な検出方法が望まれている.
我々は, 銀河の空間分布がもつ位相幾何学的情報を抽出することで BAO シグナルの検出や宇宙大規模構造を解析することを目的として, パーシステントホモロジー(PH) の計算を行った. PH はトポロジカルデータ解析(TDA) の中心的方法論であり, ホモロジー群に対して関心のあるパラメータに関するフィルトレーションによって定義される. さらに,穴の位置や形状などの宇宙論的に重要な情報を得ることができる点が特徴である. 本研究では, 共動距離(宇宙膨張の影響を取り除いた本質的な距離)をパラメータにとることにより, 空間に対する解像度を変化させた際に生じる構造を次元ごとに解析することができる.
我々はまず TDA の BAO 検出
に対する有用性の検証のため, バ
リオン有り・無しの二通りの N 体
シミュレーションデータに対してPH の計算を行った. その結果, バ
リオン有の場合にのみ有意な 2 次ホモロジー(H2), すなわち中空のシェル構造を検出した(図 3). また,
図 3: 宇宙論的シミュレーションデータから求めたパーシステン 検出された H 構造の半径の平均2ト図. 左がバリオン有, 右がなしの場合を示す.

値 は 𝑟𝑟̅ = 150.16± 8.46 [Mpc] で あり, BAO の理論と整合的であった. さらに, SDSS eBOSS サーベイによって得られた実際のクエーサー観測データに対しても同様の 解 析 を 行 い , H2 に 対 し て 𝑟𝑟̅ =146.6± 2.0 [Mpc]という推定値を得た.
また本研究では PH の逆解析も行い,得られた H2 に対応する実際の宇宙の構造を可視化することに成功した(図4). 銀河分布に対しては世界初の例である. これらの計算結果に H1, H0 についての考察を加味し, 銀河の空間分布 図 4: SDSS eBOSS データから求めた PH (左), およびそに関して得られた結果を報告する.

の逆解析によって得られた銀河分布の対応構造.

Sparse PCA by the noise-reduction
methodology
Kazuyoshi Yata
Institute of Mathematics, University of Tsukuba
and
Makoto Aoshima
Institute of Mathematics, University of Tsukuba

High-dimension, low-sample-size (HDLSS) data situations occur in many
areas of modern science such as genomics, medical imaging, text recognition,ﬁnance, and so on. In recent years, substantial work has been done on HDLSS
asymptotic theory, where the sample size n is ﬁxed or n/d → 0 as the data dimension d → ∞. Ahn et al. (2007), Aoshima and Yata (2011), Hall, Marron and
Neeman (2005), Jung and Marron (2009), Yata and Aoshima (2012) explored
several types of geometric representations on HDLSS data. Jung and Marron (2009) investigated into consistency properties both of the eigenvalues andprincipal component (PC) directions for the sample covariance matrix. Yataand Aoshima (2012) gave consistent estimators both for the eigenvalues and
PC directions together with the PC scores by developing a new PCA method
called the noise-reduction methodology. On the hand, sparse PCA (SPCA)
methods were studied by several references. Zou and Hastie (2006) and Lee,
Huang and Hu (2010) considered a regularized SPCA (RSPCA) based on L1
penalties in high-dimensional settings. Johnstone and Lu (2009) proposed athresholded SPCA (TSPCA) and showed a consistency property of the TSPCA
when n/d → 0. Shen, Shen and Marron (2013) showed that the PC directions
by the RSPCA and TSPCA have the consistency property when d → ∞ while
n is ﬁxed. In addition, Paul and Johnstone (2007) developed an augmented
SPCA method and Ma (2013) proposed an iterative thresholding procedure for
PC directions.
In this talk, we investigated the TSPCA in high-dimensional settings. We
ﬁrst proposed a new sparse estimator of the PC directions and showed that it
1holds the consistency property without any threshold values. We also proposeda shrinkage PC directions and apply it to the clustering.

References
Ahn, J., Marron, J.S., Muller, K.M., Chi, Y.-Y. (2007). The high-dimension,
low-sample-size geometric representation holds under mild conditions.
Biometrika 94, 760-766.
Aoshima, M., Yata, K. (2011). Two-stage procedures for high-dimensional data.
Sequential Analysis (Editor’s special invited paper) 30, 356-399.
Hall, P., Marron, J.S., Neeman, A. (2005). Geometric representation of highdimension, low sample size data. J. R. Statist. Soc. Ser. B 67, 427-444.
Johnstone, I., Lu, Y. (2009). On consistency and sparsity for principal components analysis in high dimensions, J. Amer. Statist. Assoc. 104, 682-693.
Jung, S. and Marron, J. S. (2009). PCA consistency in high dimension, low
sample size context. Ann. Statist. 37, 4104-4130.
Lee, S., Huang, J. Z., Hu, J. (2010). Sparse logistic principal components analysis for binary data. Ann. Appl. Statist. 4, 1579-1601.
Ma, Z. (2013). Sparse principal component analysis and iterative thresholding.
Ann. Statist. 41, 772-801.
Paul, D., Johnstone, I. (2007). Augmented sparse principal component analysis
for high dimensional data. TechnicalReport, UC Davis.
Shen, D., Shen, H., Marron, J.S. (2013). Consistency of sparse PCA in highdimension, low sample size contexts. J. Multivariate Anal. 115, 317-333.
Yata, K., Aoshima, M. (2012). Eﬀective PCA for high-dimension, low-samplesize data with noise reduction via geometric representations. J. Multivariate
Anal. 105, 193-215.
Zou, H., Hastie, R. (2006). Sparse principal component analysis. J. Comput.
Graph. Statist. 15, 265-286.

2

Decimating Nonlinear Response inState-Space Models for Battery
Capacity Estimation
Xun Shen (Tokyo University of Agriculture and Technology)
Jiancang Zhuang (Institute of Statistical Mathematics)

1

Introduction
A discrete state space model is used to perform the State of Charge (SoC) estimation, whichis written as:(1)[

Ue,t+1
SoCt+1

(2)

]
=

[
∆τ
e− Re Ce
0

]
[ ]
][
] [
∆τ
v1,t
Ue,t
0
1 − e− Re Ce
it +
+
,
ηc ∆τ
v2,t
1 SoCtCn

Ut = UOCV,0 + K0 SoCt + s(SoCt ) + it Ro + Ue,t + wt .

Here, SoCt denotes the SoC at time step t. Ut denote the terminal voltage. Ue is the voltageof the equivalent RC circuit. Re and Ce are the equivalent resistance and capacity for thecharacteristics include polarization, diﬀusion, and so on. Ro is the internal ohmic resistance. iis the charging current. UOCV,0 + K0 SoCt + s(SoCt ) is the Open Circuit Voltage of the battery.
s(SoCt ) is the unknown nonlinear part. Besides, vt = [v1,t , v2,t ]T and wt are the process andthe measurement noises which are Gaussian white noise with covariance Q and R, respectively.
In this work, we approximate s(SoCt ) by a non-parametric model and decimate it from themodel. Then, we can use the linear Kalman ﬁlter to achieve SoC estimation.

2

Proposed Algorithm
The state-space model for SoC estimation can be generally written as(3)

xt+1 = F xt + Hut + vt ,(4)

yt = C + Gxt + Jut + s(xt ).

We introduce non-parametric model to approximate s(xt ). Since s(xt ) is from Rk → Rd , we cansplit s(xt ) as
(5)s(xt ) = [s1 (xt ), ..., sd (xt )]T .

The non-parametric model used to approximate s(xt ) is deﬁned as ŝnp (xt ). We can also split(6)ŝnp (xt ) = [ŝnp,1 (xt ), ..., ŝnp,d (xt )]T .

Each ŝnp,j (xt ), j ∈ 1, ..., d can be written as(7)ŝnp,j (xt ) =Nnp∑

βij h(xt , aij , bij )i=1

where Nnp denotes the number of hidden notes of the non-parameter model, h(·) denotes theactivation function, and βij denotes the weight connecting the i−th hidden node and the output,
aij = [aij,1 , ..., aij,k ] represents the weight vector towards xt , and bij is the scalar threshold of
the i−th hidden node. βj = [β1j , .., βNnp j ]T denotes the weight vector connecting the hidden
nodes and the ouput. According to the universal approximation theorem (??), ∀ϵS > 0, ∃Npn ∈N+ , β, b ∈ RNnp , a ∈ RNnp ×k , and ∃g : Rk × RNnp × RNnp ×k → R, such that∥ŝnp,j − sj (xt )∥ ≤ ϵs , ∀xt ∈ Rk .

(8)

Then, we apply iterative algorithms to calibrate the parameter and estimate the SoC which issummarized as
(n

(n

)

)

1. Initialize iteration number nite = 0, parameter decision θ̂m ite , θ̂s ite = 0, and relaxation
factor α(nite ) = 1;
(n

)

(n

)

(n

(n

)

)

ite
as
2. Calcualte ŝnpite (x̂t|tite ) based on θ̂s ite and obtain linear part of output ylinear,t

(n

)

ite
ite )
ylinear,t= yt − α(nite ) ŝ(n(x̂t|t ) − Ĉ (nite ) − Jˆ(nite ) ut .
np

(9)

(n

)

(n

)

(For nite = 0, use ŝnpite (x̂t|tite ) = 0.)
(n

(n

)

)

(n

)

(n

)

(n

)

ite
ite
ite
ite
= {ylinear,0, ylinear,1, ..., ylinear,T}
3. Implement Kalman ﬁlter (?) based on θ̂m ite , UT = {u0 , u1 , .., uT }, Ylinear,T

(n

)

(n

)

and obtain the update of state estimation x̂t|tite , nonlinear residual ênl,tite = yt − Ĉ (nite ) −(n )Jˆ(nite ) ut − Ĝ(nite ) x̂ ite for all t;
t|t

(n

+1)

θ̂s

t

by
4. Update the parameter decision θ̂s ite∑ (n )∥ênl,tite − α(nite ) ŝnp (x̂t|t , θ̂s )∥2 .
(10)θ̂s(nite +1) = min(n

5. Update θ̂m ite
(11)

+1)

, α(nite ) by(nite +1){θ̂m
, α(nite +1) } = max log py0 (y0 ) +θ̂m ,αT∑log pyt (yt |θ̂m , α, θ̂s(nite +1) , Yt−1 , Ut−1 )t=1

or equivalently(12)(nite +1){θ̂m
, α(nite +1) } = minT∑θ̂m ,α t=1∥yt − ŷt (θ̂m , α, θ̂s∗ , Yt−1 , Ut−1 )∥2 .

∗6. Set nite = nite + 1. If nite > Nite , stop and output θ̂m= θ̂m ite , θ̂s∗ = θ̂s ite , α∗ = α(nite ) .
(n

)

(n

)

Otherwise, go back to step 2.
Gradient descent algorithm can be used to solve Eq. (10). Besides, the EM algorithm iseﬃcient to solve (12).

二部グラフ構造を持つ処理集合の実験計画とその深層学習への応用
地嵜 頌子 (大阪工業大)
藤原 良 (筑波大学)
宮本 暢子 (東京理科大学)
本報告では, 完全二部グラフの辺を処理集合に持つようなブロック計画を提案する. 実験計画では普通は処理の集合とブロック (処理の部分集合) の集まりでブロック計画を構成する. ここでは処理集合は単なる集合ではなく, その中に構造を持っている場合を考える. いま完全二部グラフの辺の集合が処理集合であるとする. そのとき処理パラメータの推定精度を上げるために処理のペアをバランスさせようとすると, 複雑な問題が生じる. 完全二部グラフ Kn,m の二つの頂点集合を

V1 , V2 とし, 以下のような条件を満たす部分グラフの集合 B = {B1 , B2 , . . . , BN } を考える.
条件(1) B のどの部分グラフ Bi も V1 , V2 の点を全て含む(2) B の中に Kn,m のすべての辺はちょうど µ 回ずつ現れる(3) 任意の x1 , x2 ∈ V1 , (x1 ̸= x2 ), y1 ∈ V2 に対して, 二つの辺 (x1 , y1 ), (x2 , y1 ) が同時に含まれる部分グラフは B 内に λ21 個存在する(4) 任意の x1 ∈ V1 , y1 , y2 ∈ V2 , (y1 ̸= y2 ) に対して, 二つの辺 (x1 , y2 ), (x1 , y2 ) が同時に含まれる部分グラフは B 内に λ12 個存在する(5) 任意の x1 , x2 ∈ V1 , (x1 ̸= x2 ) かつ y1 , y2 ∈ V2 , (y1 ̸= y2 ) に対して, 二つの辺 (x1 , y1 ), (x2 , y2 )が同時に含まれる部分グラフは B 内に λ22 個存在する
これらの条件を確認するために部分グラフ B1 , B2 , . . . , BN に対応する (n×m)-結合行列をそれぞれ H1 , H2 , . . . , HN とする. 結合行列 H1 , H2 , . . . , HN に対して, (N × m)-行列 Xj (j = 1, 2, . . . , n)
をその第 i 行が Hi の第 j 行となっている行列であるとする.

それらを結合した (N × nm)-行列

X = [X1 |X2 | · · · |Xn ] を計画行列という. このとき, 上の条件 (2) から (5) が満たされれば, 情報行
列 X t X の小行列に次のように現れる: (a) 小行列 Xit Xi の対角成分は全て µ (条件 2), 非対角成分は全て λ21 (条件 3). (b) 小行列 Xit Xj , (i ̸= j) の対角成分は全て λ12 (条件 4), 非対角成分は全て λ22 (条件 5). 理想的には λ21 = λ12 = λ22 (= ρ), すなわち X t X = µI + ρ(J − I) の形になる
時, 推定量の精度が最適になる. そうでない場合には, 非対角成分に現れる値の種類ができるだけ少ない方が良いとされている. (長尾, 栗木 [3])定義 1 ((r, λ)-design [2]) V を m 点集合とし, B を V の部分集合からなる集合族 (ブロック集合と呼ぶ) とする. (V, B) が次の条件を満たすとき, (r, λ)-design という.

• V の各要素は必ず B の r 個のブロックに含まれる• V の任意の異なる 2 つの要素を同時に含むブロックは, B の中に必ず λ 個存在する(r, λ)-design は点の数を m, ブロック数を n とすると, (n × m)-結合行列 H で表現できる. そして,
V 上の巡回写像を P0 , P2 , . . . , Pm−1 であるとし, m 個の design としては同型な結合行列を生成することができるとする. これらを基本結合行列とする.

Hi = HPi , i = 0, 1, . . . , m − 1
これらの基本結合行列からできる (m × nm)-計画行列 X = [X0 |X1 | · · · |Xn−1 ] は条件の (1) から(4) を満たす. ただし (5) の条件は満たさない. (5) の条件を満たすために, 次の概念を導入する.
1

定義 2 (Diﬀerence matrix [1]) Zm 上の (k × mγ)-行列を D = [dij ] とする. 任意の 1 ≤ i < j ≤ kについて, Zm 上の多重集合{dih − djh(mod m) | h = 1, 2, . . . , mγ}

が Zm の要素をちょうど γ 回ずつ含むとき, D を (m, k, γ)-diﬀerence matrix (DM) という.
先ほど巡回的に作られた Hi = HPi (i = 0, 1, . . . , m − 1) を基に計画行列 X = [X0 |X1 | · · · |Xn−1 ] を DM を用いて次のように拡大する. ここで, Xi は各 Hj (j = 0, 1, . . . , m − 1) の i 行目からなる行列である. (Hi , Xi の行列の行を, 上から 0 行目, 1 行目, ... と数える.)



Xd11







 Xd 
12
X1 = 
,
..


.


Xd1 (mγ)

Xd21





 Xd 
22
X2 = 
,...,..
.

Xd2 (mγ)



Xdk 1





 Xd

k
2


Xk = 
.
..


.


Xdk (mγ)

ただし, D = [dij ] は (n, k, γ)-DM とする. このとき, 拡大された計画行列 X = [X1 |X2 | · · · | Xk ]は条件 (2) から (5) を満たす. 拡大された計画行列 X でも, 各行は一つの実験パターンを表している. 拡大前の計画行列 X の第 i 行は Hi の結合行列に対応していたが, 拡大された場合は変形された結合行列の集合になっている. その時の新しい結合行列のサイズは k × m となっている, (これを実験結合行列と呼ぶこととする).
ここでパラメータを整理すると,• 基本結合行列 Hi のサイズは n × m• Hi を生成する巡回群の位数は m• Diﬀerence matrix (m, k, γ)-DM のサイズは Zm 上の k × mγ
• 拡大計画行列のサイズは m2 γ × mk, この時の実験結合行列のサイズは k × m• (r, λ)-design と (m, k, γ)-DM を使って情報行列を作った場合, µ = rmγ, λ21 = mλγ, λ12 =λ22 = µγ + λγ(m − 1) なる計画行列を作ることが可能である.
本報告では, 応用事例や具体例を示すとともに, diﬀerence matrix の要素や基本結合行列を生成する際に用いる巡回群について議論した.

謝辞

本研究は JSPS 科研費 19K11866 の助成を受けたものです.

参考文献[1] T. Beth, D. Jungnickel, and H. Lenz. Design Theory. Encyclopedia of Mathematics
and its Applications. Cambridge University Press, 2 edition, 1999.
[2] D. Stinson. Combinatorial designs: constructions and analysis. Springer Science &Business Media, 2007.
[3] 長尾壽夫 and 栗木進二. 数理統計学. データサイエンス・シリーズ. 共立出版, 2006.

2

マクロ進化の⼤規模種間⽐較研究における枝特異的⽅向性淘汰モデルの提案

⼤久保祐作 12
1. 統計数理研究所モデリング研究系 特任研究員
2. データサイエンス共同利⽤基盤施設データ同化研究⽀援センター 特任研究員
・背景

種間系統⽐較法(phylogenetic comparative method; PCM)とは⽣物形質のマクロ進化を分析す
る統計的⼿法の総称である。PCM は進化⽣物学で提案され(Felsenstein 1985)、祖先種の状態や進化速度などのパラメータ推定や、回帰分析において進化史の共有に起因する相関を補正したりする⽬的で広く応⽤されている。特に後者については、環境変動の与える⽣物への影響を評価する際に複数の種を跨ぐ⼤規模なデータセット(例えば > 1000 種)でメタ解析を⾏う⼿法が広く普及したことに伴い近年急速に利⽤が広がっている(中川 & 久保 2016)。
既存の PCM には⼤きく分けガウス過程によるアプローチとシミュレーションによるアプローチの 2 つが知られてきた。ガウス過程によるアプローチでは、過去から現在に⾄る⽣物の進化を Wiener 過程や Ornstein–Uhlenbeck 過程などのよく知られた確率過程としてモデル化する。このアプローチは推定量の信頼性が⽐較的評価しやすく、計算も容易なことから⻑く
PCM で主流であった。しかし、モデルの柔軟性が不⼗分で実データへの適合や予測精度に問題があることが指摘されている(Losos 2011)。また、モデルの構造にミクロ進化など他の⽣物学領域で知られた知⾒との乖離があり推定結果の解釈が困難であることも知られてきた(Pennell et al. 2013)。特にデータセットが⼤規模になると多様な進化パターンを遂げた種を内包するため、⽣物学的な背景理論を取り⼊れることのできる既存のモデルより柔軟性が⾼いモデルが求められてきた。
⼀⽅シミュレーションによるアプローチでは、集団遺伝学など関連領域の知⾒を取り⼊れた柔軟なモデルが可能となるが、閉じた形式の尤度が得られないことが多い。例えば
Kutsukake & Innan(2013; 2014)は系統樹上の⼀部に特異的な⽅向性選択が⽣じる過程をシミュレーションし、近似ベイズ計算(Approximate Bayesian Computation; Beaumont et al. 2002)を⽤いて計算する⼿法を提案した。この⼿法はガウス過程アプローチに⽐べ極めて柔軟な⼿法として評価されてきた(Pennell et al. 2013; Garamszegi2014 など)⼀⽅で、尤度ベースの⼿法と⽐べ 1) 統計的信頼性が評価しにくい 2) ⾼い計算コストを伴う 3) 現状では各データセ
ットに合わせた専⽤ソフトウェアが必要で、実務家にとって利⽤しづらい等の問題がある。
・提案⼿法

そこで本研究では、Kutsukake & Innan(2013; 2014)のシミュレーションモデルをガウス過程
モデルとして再定式化した“枝特異的⽅向性淘汰モデル(branch-specific directional selection;

BSDS)モデル”を提案する。このモデルは 2 次オーダーで元のシミュレーションと等価で、⾼次は進化史が⼗分⻑ければ 0 に収束することを明らかにした。
・数値実験

ランダムに⽣成された系統樹のもとで Kutsukake & Innan(2013; 2014)のモデルに従う⽣物進化のシミュレーションを⾏い、⽣成されたデータに提案⼿法を適⽤した。その結果、祖先種における形質値、進化速度、淘汰の強さなどのパラメータを不偏に推定できることを明らかにした。⼀⽅、既存のガウス過程モデルと近似ベイズ 計算による⽅法では、それぞれ共通
祖先形質と淘汰圧を過⼤に推定することを明らかにした。
・実データ分析

霊⻑類のデータを⽤いて、ヒト(Homo sapience) の脳サイズがチンパンジー(Pan troglodytes)と分岐して以降にどの程度の淘汰圧を受けたか推定した。その結果、脳サイズが⼤きくなる
⽅向に⽣じた変異が他の種に⽐べ 1.73 倍 (95%確信区間 [1.62-1.87])多く⽣じており強い⽅向性をもった淘汰が⽣じたと推定された。
・今後の課題本研究では系統樹上の⽅向性淘汰を⽣じた枝が既知であるという仮定で淘汰圧を推定する⼿法を検討し、提案した BSDS モデルが不偏な推定を可能にすることを明らかにした。しかし、より⼤規模なモデルについては事前に変数選択が必要になるなど、この仮定が不適当な場合がありうると考えられる。
⼀般に⼤規模なモデルを扱う場合には、近年スパース推定により必要な変数の選択とパラメータ推定を同時に⾏うアプローチが有効であると知られている。ところが PCM の場合、系統樹という複雑に構造化されたモデルを扱うため、よく知られた LASSO や Group LASSOなどの⼿法では良い推定を与えない。今後、系統樹の構造を適切に反映した罰則項の検討や、その罰則項に対応する事前分布のモデル化などを検討する必要がある。

2 次林再生における林冠個体の 15 年間変遷を観る統計モデル
島谷健一郎（統計数理研究所）
本研究では、1970－80年代に青森に設定されたブナ天然更新指標林で得ている、2005
年からの15年間の樹木の動態データを分析した。
幅5m, 長さ170mの調査地を設け、その中の胸高直径2.5cm以上の高木種について、樹種、胸高直径、樹高、階層（目視により１＝林冠層、２＝亜林冠層、３＝中―低木層、
４＝倒れなど、５＝枯死）を4回、計測した。樹木の位置は、170mを5mごとに区切ったコドラート番号で与えた。計測したのは約900本である。その樹高と直径の変化を以下に示す。
樹高分布の変化

直径分布の変化

250

2-3cm
6
10
14

250

2004
2013
200

2017

200

3-4cm
7
11
15

4
8
12
16

5
9
13
17

個体数

個体数

150

100

150

100

50

50
0
100

200

300

400

500

600

700

樹高 (cm)800

900 1000 1100 1200

0
2004

2013

2020

概して針葉樹、特に人工林では、最初に高かった個体が良好な光環境により成長が促
進され、最初に出遅れた個体はその後塵を拝する。しかし、天然性広葉樹林では、その限りでなく、以下のような逆転現象も観測される。
個体 2 と個体 3 は同じくらいの樹高だったが、直径で劣っていた個体２が樹高でも直径でも階層でも先を進んでいる。個体 1 は個体２より高く直径はほぼ同じだったが、樹高でも階層でも後塵を拝している。

当初多数観察された林冠層の個体は、年と共に競合により亜林冠層あるいは枯死といった運命をたどった。

2004 年の林冠木の 2020 年の運命

階層分布の変化

個体数

350

100%

300

90%

250

80%

200

70%

枯死

60%

倒れ

50%

低⽊層

40%

亜林冠

30%

林冠⽊

150
100
50
0
2005

Jul-05

2013

2016

2017

2020

20%
10%

林冠木

亜林冠層

低木層

倒れ

枯死

0%

本研究では、直径－樹高の関係式（アロメトリー）に着目した。ずっと林冠層だった
ブナ111本を抽出した。直径と樹高については古くから様々な回帰式が提案されている。
その主な目的は材積、バイオマス、2酸化炭素吸収量の推定などに置かれている。そのため、樹高測量はスナップショットで広い森林の代表的な樹高の個体に絞って行われることも多い。一方、本調査では、ブナ2次林の更新過程を詳細に追跡するべく、個体単位で樹高を正確に測量してきた。
直径 (x)－樹高 (y)アロメトリーでしばしば使われる数式に、Ogawa (1965) がある。
1
1
1


y ax b H Max

Hmaxは潜在的にとりうる樹木の樹高の最大値である。
この式は、単に当てはまりの良さから使われているのではなく、この式が、直径と樹高それぞれの相対成長が比例する、ただし樹高の成長は最大値が近づくにつれて漸減す
る、という以下の微分方程式から導かれるという特徴がある。

1 dy
1 dx
y
b
(1 
)
y dt
x dt
H Max
そこで、本研究では、5回の計測データを用いて、結果としての回帰式に加え、相対成長の関係についても吟味した。
結果として、相対成長比例定数 (b) は個体及び年による変動が見られ、上式が成り立っているとは言い難いが、決して見当違いな関係式にも見えなかった。個体や環境（金林個体との競合、密度など）による変異を適切に分散パラメータとして表現する階層モデル（Ogawa (1965) 式を状態方程式とし樹高などを観察モデルとする状態空間モデル）で定式化する方向性が伺えた。