科研費（基盤A 20H00576)「大規模複雑データの理論と方法論の革新的展開」
研究代表者：青嶋 誠（筑波大学）

予測モデリングの理論と応用

開催責任者：小森理（成蹊大学）・深谷

肇一（国立環境研究所）

日時：2021 年 10 月 22 日（金）～24 日（日）場所：成蹊大学 6 号館 401 室（ハイブリッドで実施）
内容・目的：予測モデリングは理論と応用の両面で非常に重要であり，医学，生態学，生物学，工学，ビジネスなどの分野で盛んに研究されているテーマの一つです．特に近年注目されている大規模かつ複雑データを扱う予測モデリングの研究は未開拓な部分が多く，これから大いに進展する分野であることが期待されます．本シンポジウムでは様々な分野で活躍される研究者の講演を広く募集し，分野横断的な議論を通し新たな研究の方向性を探っていければと思います．
旅費の配分：講演者を中心に配分します。旅費希望の方は講演申込の際にお伝え下さい。
宿舎の斡旋：斡旋しません。
講演申込期限：2021 年 9 月 24 日（金）氏名・所属・講演題目を電子メールでお知らせ下さい。
予稿期限：2021 年 10 月 8 日（金）A4 サイズ 10 頁以内で作成し、PDF ファイルを送信して下さい。
報告書原稿：
報告書を作成しますので、予稿とは別に報告書原稿（A4 サイズ 2 枚）も PDF ファイルで送信して下さい。
問い合せ先・講演申込先・予稿送付先・報告書原稿送付先：〒180-8633 東京都武蔵野市吉祥寺北町 3-3-1
成蹊大学

理工学部情報科学科

小森理

Email: komori[at]st.seikei.ac.jp TEL: 0422-37-3764

プログラム10 月 22 日（金）特別講演
13:00-13:50 岡村寛(中央水産研究所)

水産資源データと予測モデル
13:50-14:25 Dou Xiaoling(早稲田大学データ科学センター)
Estimation of B-spline copula
14:25-15:00 米岡

大輔(聖路加国際大学)

予測モデルのメタアナリシス
15:00-15:20 休憩
15:20-15:55 大前勝弘(国立循環器病研究センター)

マジョリティーに対する治療効果のロバストな予測
15:55-16:30 三角俊裕*（横浜市立大学）, 三枝

浜市立大学), 藤田豊（横浜市消防局), 古場祐輔（横浜市立大学), 原亮介(横浜市消防局）

横浜市における救急出場件数に対する関数動的予測モデリング
16:30-17:05 得丸

久文(カラハリプロジェクト)

概念の複雑性・信頼性予測のための受け入れ検査(デジタル言語学)

みゆひ（横10 月 23 日（土）午前
特別講演
10:00-10:50 江口

真透(統計数理研究所)

予測のための最大エントロピーと最小ダイバージェンス
10:50-11:25 野津昭文(静岡県立静岡がんセンター)

臨床研究における統計家の役割
11:25-12:00 岡山

大成*(千葉大学融合理工学府), 内藤

貫太(千葉大学理学研究院)

ガンマクレームによる資産過程における破産確率の漸近近似
12:00-13:00 お昼休み
特別講演
13:00-13:50 久保田康裕(琉球大学，株式会社シンクネイチャー)

生物多様性ビッグデータと地球環境保全
13:50-14:25 楠本聞太郎(九州大学 農学研究院)

生物分布の時空間予測：機構論的モデル構築に向けた課題
14:25-15:00 城田慎一郎(明治大学商学部)

Preferential Sampling に関する最近の話題について
15:00-15:20

休憩

15-20-15:55 矢島

豪太*(日本大学生物資源科学部), 中島 啓裕(日本大学生物資源科学部)

近似ベイズ計算による野生動物の集団サイズ推定 - 生成過程が複雑なデータにどう対処するか？ 15:55-16:30 深谷

肇一(国立環境研究所)

環境 DNA による生物多様性評価のための階層モデルとベイズ決定分析
特別講演
16:30-17:20 島津

秀康(Loughborough University)

生物多様性変化の定量化とその理解

10 月 24 日（日）特別講演
9:30-10:20 島谷健一郎(統計数理研究所)

定量的生態学と統計モデル：食われる側からみた特定食う種の貢献度
10:20-10:55 川森

愛(統計数理研究所)

ニワトリ雛の採餌行動解析ー理論最適解から乖離する理由を探すー
10:55-11:30 永井勇(中京大学 教養教育研究院)

ランク落ちでの無罰則推定法と不偏推定量
11:30-11:55 小森理(成蹊大学)A unified formulation of k-means, fuzzy
Kolmogorov-Nagumo averagec-means and Gaussian mixture model by the

⽔産資源データと予測モデル
⽔産研究教育機構

⽔産資源研究所
岡村

漁業情報解析部

寛

海洋⽣物の再⽣産関係（親⼦関係）データは，不確実性が⼤きいことで知られ，再⽣産関係の推定には頑健推定と呼ばれる⽅法が使⽤されることがある（Chen and Paloheimo 1995,
Shertzer and Prager 2002）．我が国の⽔産資源でもいくつかの種（系群）に対して頑健推定法が使⽤されている．⼀⽅で，再⽣産関係データは時系列データであり，⾃⼰相関の推定が重要である．頑健推定は外れ値のような“良くない”データに「鈍感」である必要があり，⾃⼰相関推定はデータ間の関係性に基づくので他のデータに「敏感」である必要がある．つまり，
2 つの⼀⾒⽭盾した問題への対応が必要となるが，その両⽅をうまく取り扱う頑健推定法の開発はこれまで⼗分になされてこなかった．そこで，頑健推定法を拡張し，外れ値に強く，かつ⾃⼰相関の正確な推定も可能となる推定法を考案した（Okamura et al. 2021）．
加⼊尾数 R の対数値を r = log(R)と書き，親⿂量を S で書くとすると，加⼊尾数は再⽣産
関係 f(S|θ)の周りの誤差 ε をあわせて𝑟𝑓 𝑆 |𝜃

𝜀で与えられる（t = 1, …, T）．このとき，残差 εt は，
𝜀 ~𝑁 𝜌 𝜆𝜀, 𝜎 ⁄𝜆という形であるとする．ここで，λt を，
𝜆

exp

𝜙𝜀

と定義する．1–λt は，各データが外れ値であるかどうかの確率に対応し，上の⾃⼰相

図 1. 時系列クロスバリデーションの模式図．1

関の式から，λt が⼩さくなると（外れ値で年ずつデータを除いたセット（青色部分）で予ある確率が⼤きくなると）⾃⼰相関に対す測モデルを構築し，除かれたデータ（白色部る寄与が⼩さくなることになる．逆に，λt が分）を予測する

⼤きい場合は，⾃⼰相関への寄与が⼤きくなる．𝜙 → 0とすれば，通常の AR(1)モデルとなる．また，この確率分布は，ρ=0 とすると，Tukey の biweight 関数と似たような形になる．
𝜙は外れ値の程度を決めるパラメータであり，各データが再⽣産関係推定に⼤きく寄与するものか，外れ値である可能性が⾼いのかを決定するための鍵となるものである．そこで，𝜙パラメータを加⼊尾数に対する時系列クロスバリデーションによって決定してやることにする（図 1）．𝜙以外のパラメータを最尤推定によって求めたとき，最もクロスバリデーションエラーが⼩さい𝜙の結果を選択する．

図 2. シミュレーション結果の一例．大きな外れ値があり，自己相関も高い（ρ=0.8）の場
合．赤いボックスは長期加入予測，緑のボックスは短期の加入予測に対応する．最小絶対値法（LAD）と最小二乗法（LS）はバイアスした結果となるが，新しい頑健推定法（RSR）は不偏に近い推定ができている
⼤きな外れ値や⾃⼰相関を想定したシミュレーションによって，加⼊量の短期・⻑期両⽅の予測の観点で新しい⽅法は従来の⽅法より良いパフォーマンスを持つことが⽰され（図2），実データに適⽤した結果も変わる可能性があることが⽰唆された（図 3）．理解志向型のモデルと応⽤志向型の統計⼿法（江崎 2020）を効果的に組み合わせることが，⽔産資源
データの解析と理解には不可⽋なことである．頑健推定と⾃⼰相関推定の問題は，再⽣産関係推定にとどまるものではなく様々なデータ解析に現れるものであることから，新⼿法は広い応⽤範囲を持つものと考えられる．
引用文献
Chen, Y. and Paloheimo, J. E. 1995. A
robustregression

analysisofrecruitment in fisheries. Can. J.
Fish. Aquat. Sci. 52, 993–1006.
江崎

貴裕 2020. データ分析のための数理モデル⼊⾨ 本質をとらえた分析のために．ソシム．
Okamura, H., Osada, Y., Nishijima, S.,
and Eguchi, S. 2021. Novel robust
time series analysis for long-term

図 3. マアジ太平洋系群に新しい頑健推定モデルをand適用した結果short-term

prediction.

Scientific Reports 11: 11938.
Shertzer, K. W. and Prager, M. H. 2002. Least median of squares: a suitable objective function forstock assessment models? Can. J. Fish. Aquat. Sci. 59, 1474–1481.

Estimation of B-spline copula
Xiaoling Dou*, Satoshi Kuriki† , Gwo Dong Lin‡ , Donald Richards§

B-spline copula is constructed with B-spline basis functions. It is a generalization of the Bernstein
copula, and its form is similar to that of the Bernstein copula (Dou, et al., 2021). For given p + 2d + 1
knots
s−d = · · · = s0 = 0 ≤ s1 ≤ · · · ≤ sp−1 ≤ 1 = sp = · · · sp+d ,we can generate d + p (=: n) B-spline functions on [0, 1] of degree d,i = −d, . . . , −1, 0, 1, . . . , p − 1.

Nid (s),For any u, v ∈ [0, 1], the bivariate B-spline copula density is deﬁned asc(u, v; R) =

m ∑
n
∑

rkl ϕk,m (u)ϕl,n (v)(1)

k=1 l=1

whererkl ≥ 0,n∑rkl = qk,m ,

m
∑

and
ϕk,m (s) =rkl = ql,n ,

k=1

l=1

∫ 1

1
qk,m

d
Nk−d−1(s),dNk−d−1(s)ds.

qk,m =
0

When m = n, it is known that the maximum correlation of the B-spline copula c(u, v; R) is attainedwhen
R = (rkl ) = diag(qk,n )1≤k≤n ,and we denote the B-spline copula with the maximum correlation asc+n (u, v) =n∑

qk,n ϕk,n (u)ϕk,n (v),u, v ∈ [0, 1].

(2)

k=1

Suppose that we have a set of i.i.d. data (xt , yt ), t = 1, . . . , N . We consider methods of estimating
the joint density by using the B-spline copula. As the B-spline copula has a similar form to that of theBernstein copula, the EM algorithm methods for the Bernstein copula (Dou, et al, 2016) are applicable
after some changes made on the basis functions and the constrictions on the parameter matrix R. We
consider two EM algorithms for the two forms of the B-spline copula in (1) and (2), respectively. In bothmethods, the marginal distributions and the copulas are estimated separately. That is, for each marginal,
we use the kernel method and the empirical distribution function to estimate the density function and thecumulative distribution function, respectively.

* Waseda University, † The Institute of Statistical Mathematics,‡ Academia Sinica, § Pennsylvania State University

1

To estimate the joint density of two random variables, the general idea is to use the general form (1)directly,h(x, y; R) = c(F (x), G(y); R) f (x)g(y).
In this method, we need to determine the parameter matrix R and its size m × n. The ﬁrst EM algorithm
is to estimate the parameter matrix R of the ﬁnite mixture distribution (1). The size of R, m × n, can be
determined by AIC.
When the two random variables X and Y are highly correlated, we can use the copula density of the
maximum correlation case (2) and the independent copula to approximate the joint density of the databyh(x, y; q, n) = (1 − q)f (x)g(y) + qc+
n (F (x), G(y))f (x)g(y).
In this case, we need to ﬁnd parameters q and n. For this model, except for another EM algorithm, agrid method can also accomplish this purpose.

References[1] Dou, X., Kuriki, S., Lin, G. D., and Richards, D. (2016). EM algorithms for estimating the Bernstein copula, CSDA, 93, 228–245.
[2] Dou, X., Kuriki, S., Lin, G. D., and Richards, D. (2021). Dependence properties of B-spline copulas, Sankhya A, 83, 283–311.

2

予測モデルのメタアナリシス
Meta-analysis of prediction models
聖路加国際大学 米岡大輔
概要 本発表では、多変量メタアナリシスの手法を用いた回帰モデルの係数統合とそれにまつわる問題点を整理する。 また、発展型として現在取り組んでいる予測モデルの線形予測子の統合についても時間があれば議論する。

1

はじめに

医学、 疫学、 教育学を始めとした学際的分野において、 各研究が報告する要約統計量を統合する試
みとしてメタアナリシスがある。 従来のメタアナリ
シス研究は、 例えば各研究の疾患割合などの単一指標の統合に取り組んできた。 しかし近年、 多変量メ
タアナリシスの手法が発展し複数の指標の統合が盛んに研究されるようになり、 ネットワークメタア
ナリシスなどへの応用も見られる。 この文脈においてBecker et al. (2007)[1]は回帰モデルの係数統合を提唱し、 予測モデルのメタアナリシス的統合という分
野に脚光が当たり始めた。
本発表では、 回帰モデルの係数統合における様々な問題点（例えば、 報告されない要約統計量の問題
や、 予測タスクは同じだが回帰関数が異なる場合など）を整理し、 特に回帰モデル間に共変量の組の相違がある場合のバイアスとその補正について議論する。 また、 時間があれば発展型として現在取り組んでいる予測モデルの線形予測子の統合についても議論する。

2

提案方法: 共変量の組に相違がある場合の係数統合

だしpXZ は共変量XとZの同時分布）
。
!#
"

exp(β0∗ + β1∗ X)
1
=0
E
Y −
1 + exp(β0∗ + β1∗ X)
X
特にXとZがともに連続変数で、Z|Xが正規分布となるときに、 このgはCramer (2003) [2]やChao et al.
(1997) [3]の以下のような結果と一致することが示せる。
α0 + α2 γ0β0∗ = g0 (α0 , α1 , α2 , pXZ ) ≈ p
1 + c2 α22 σ 2
α 1 + α 2 γ1β1∗ = g1 (α0 , α1 , α2 , pXZ ) ≈ p
1 + c2 α22 σ 2
ただし、ここではZ = γ0 +γ1 X+τ、τ |X ∼ N (0, σ 2 )と仮定している。 共変量が全て連続の場合は、 このよ
うにgを陽に書けることがあるが、 一般には難しく、
上のスコア関数の普遍性の要件から陰関数gとして定義する。
このgを用いて、以下の一般化最小二乗法を解くことで、 統合された係数であるα̂∗ = (α0∗ , α1∗ , α2∗ )を得る。
arg minα

手元にK個(k = 1, . . . , K)の過去研究があり、 その各々がロジスティック回帰の係数ベクトルβk とその分散共分散行列Σk を報告しているとする。 各モデルには異なる組の共変量が入っているため、 単純に多変量メタアナリシスで統合することは難しい状況を考える。 また、 手元には一つだけ実データもあり、過去研究のすべての共変量が含まれるとする。
簡単化のためK = 2とし、 一つが真モデル(1)を推定しているのに対し、 もう一つが誤特定のモデル(2)を使った場合の係数統合について考える。
logitP (Y = 1|X, Z) = α0 + α1 X + α2 Z

(1)

K n
o>
n
o
X
θ̂k − fk (α, p̂XZ )
Σ−1
θ̂
−
f
(α,
p̂
)
k
k
XZ
k
k=1

ただし、 p̂XZ は手元のデータより推定した量であり、α̂ = (α̂ , α̂ , α̂ )> (if k = 1)
k
0
1
2
θ̂k =
β̂k = (β̂0 , β̂1 )> (if k = 2)fˆk (α, p̂XZ ) =

3


α̂

k (if k = 1)
(g1 (α, p̂XZ ), g2 (α, p̂XZ ))> (if k = 2)

数値実験

K = 9とし、うち3研究は共変量X, Z、３研究はXのみ、 ３研究はZのみでロジスティック回帰を行った
αとβのパラメータの関係を表すβi∗ = gi (α0 , α1 , α2 , pXZ ) 結果の係数統合を考える。 比較対象として、 最初の フ ル セ ッ ト の3研 究 の み の 多 変 量 メ タ ア ナ リ シ(i = 0, 1)は以下のスコア関数の普遍性を満たす（たlogitP (Y = 1|X) = β0 + β1 X

(2)

ス（M1）、 脱落している係数を平均値補完した後の多変量メタアナリシス（M2）を用意し、 バイアスとMSEを1000回のモンテカルロシミュレーションにより評価する。 また、 簡単化のために（脱落しているものも含めて）共変量はすべて正規分布に従うとする。
結果は以下のTableに示す通りである。 M1やM2と比較して、バイアスは10-50%程度の改善、MSEは1090%程度の改善が見られた。

[3] Chao WH, Palta M, Young T. Effect of omitted
confounders on the analysis of correlated binary
data. Biometrics 1997; 53: 678-689.

MSE

バイアス

-3
提案手法
0.133
0
-0.377
1
0.149
2

4

True value of
-1
0

1

1

3

0.100
0.070
0.080

0.166
0.350
0.138

-3
提案手法
0.058
0
0.251
1
0.049
2

従来手法: フルセットのみ (M1)
0.114 -0.005 -0.022 0.012
0-0.1400.016 -0.001 -0.003
1
0.066 -0.012 -0.014 0.001
2
従来手法: 平均値補完 (M2)
0.687
0.375
0.204
0
-0.915 -0.255 -0.013
1
0.756
0.379 -0.015
2

0.062-0.060
0.070

0.048
0.000
0.036

0.378
0.223
0.326

True value of
-1
0

1

1

3

0.030
0.016
0.021

0.074
0.239
0.048

0.112
0.131
0.058

従来手法: フルセットのみ (M1)0.101 0.042 0.030 0.046
00.188 0.018 0.003 0.018
1
0.042 0.021 0.017 0.022
2

0.105
0.183
0.042

0.687
0.857
0.688

従来手法: 平均値補完 (M2)0.480 0.149 0.049 0.151
00.888 0.070 0.001 0.055
1
0.576 0.148 0.008 0.111
2

0.481
0.780
0.477

結論

本研究では、 メタアナリシスにおける回帰モデルの共変量の組の研究間相違の問題を、 脱落変数バイスの問題と捉え直し、 そのバイアスをスコア関数の普遍性の要件を考えることで解決した。 更にシミュレーションと実データを用いて、 既存手法よりバイ
アスやMSEの意味で本手法が改善されることも示した。

参考文献[1] Becker BJ, Wu M-J. The Synthesis of Regression
Slopes in Meta-Analysis. Statistical Science 2007;
22: 414-429.
[2] Cramer JS. Logit models from economics and otherfields. Cambridge University Press: Cambridge,
UK ; New York, 2003.

0.024
0.015
0.019

0.013
0.002
0.008

マジョリティーに対する治療効果のロバストな予測
国立循環器病研究センター

大前勝弘

ある疾患を有する対象集団に対して，2 つの治療選択 a ∈ A = {0, 1} のうち，どちらがより患者の状態を改善するかを予測したい．このために，被験者の治療前の状態ベクトル x ∈ Rp に応じて治療を決定する関数 d : Rp → {0, 1} を考える．この d を治療方針と呼ぶ．
患者の状態を表す確率変数を Y とする．治療前の状態ベクトル x と治療 a を条件付けた下での Yの期待値を Q(x, a) = E[Y |X = x, A = a] と書くと，d(Q(x)) = argmaxa Q(x, a) = I(Q(x, 1) >Q(x, 0)) が自然な最適治療方針となる．ここで，I は指示関数である．つまり，条件付き期待値
Q(x, a) = E[Y |X = x, A = a] に関する推論を行うことで，最適治療方針を推定することができる．
医学研究で得られる質の高いデータはそれほど膨大でないことと，治療選択問題における解釈性の重要さから，線形モデルを用いた治療方針の推論が有望である．例えばアウトカムが連続値で正規性を仮定できる場合には，治療方針の推定問題は，正規線形交互作用モデル

Y |(X = x, A = a) ∼ N(Q(x, a; β), σ 2 )Q(x, a; β) = β0 + β1⊤ x + β2 a + β3⊤ xaのうち関心のあるパラメータ (β2 , β3 ) に関する推論の問題へと帰着される．
ところが，多くの疾患と治療に関して，治療反応性は患者ごとでしばしば異なることが知られている．例えば {(x, a = 1, y), Q(x, 1) > Q(x, −1) > y} を満たす患者のように，最適治療方針による治療選択が不合理になる場合がある．治療反応の異質性を考慮せずに対象集団全体にとって良い治療方針を模索しようとすると，このようにイレギュラーな特性をもつ学習サンプル：マイノリティーによって歪められた推論になってしまう可能性がある．そのような推論をもとに治療選択をすると，大多数の標準的な特性を持つ患者集団：マジョリティーにとって不利益を生じうる可能性もある．治療開発のコストパフォーマンスを高め，さらなる治療開発へと繋げるためには，「マジョリティーに対して有効な治療方針の推論」および「推定された治療方針が有効な集団の特定」が重要な課題となる．
本発表では，これらの問題に対応するために，治療周辺回帰モデルpθ (y|x) =∑pθ (y|x, a)p(a|x)

(1)

a∈A

を提案した。この (1) 式に基づくパラメータ θ = (β0 , β1 , β2 , β3 , σ 2 ) の対数尤度関数は，ℓR (θ) =

n
∑
i=1

(
log

∑

)
p(yi |xi , a; θ)p(a|xi )a∈Aで与えられる．従って，治療周辺回帰モデルの関心のあるパラメータに関する推定方程式は，次式のような正規線形回帰モデルの推定方程式 EF (β, x, a, y) =y−Qβ (x,a) ∂σ2∂β Qβ (x, a) の事後確率重み付き

図 1 正規線形回帰モデルと治療周辺回帰モデル

平均となる．

EM (β, x, a, y) = pβ (−1|x, y)EF (β, x, −1, y) + pβ (1|x, y)EF (β, x, 1, y).

(2)

ただし，





1
1 + exp(−α1 (x)y − α0 (x))
pβ (a|x, y) =1

1 + exp(α1 (x)y + α0 (x))

if a = 1
if a = −1であり，α1 (x) =Qβ (x, 1) − Qβ (x, −1)Qβ (x, −1)2 − Qβ (x, 1)2
p(1|x)
,
α
(x)
=
+ log
.
0
2
2
σ
2σ
p(−1|x)正規線形回帰モデルの最尤推定量はマイノリティー集団に敏感であり，マイノリティー集団の存在によってパラメータの推定が非常に不安定になってしまう．これに対して，周辺回帰モデルの最尤推定
量は，治療レジメンに不適合なマイノリティー集団の存在に対して頑健な性質を持つ（図１）．これは推定方程式の以下の性質により得られる．

suppβ (1|x, y)∥EF (β, x, 1, y)∥ < ∞y≤Q(x,0)lim pβ (1|x, y)∥EF (β, x, 1, y)∥ = 0

y→−∞

sup

pβ (0|x, y)∥EF (β, x, 0, y)∥ < ∞y≥Q(x,1)

lim pβ (0|x, y)∥EF (β, x, 0, y)∥ = 0.

y→∞

治療周辺回帰モデルは一般化線形モデルと同様の拡張が可能であるため，様々な種類のアウトカムを改善するための治療レジメン推定問題への応用が今後期待される．

2

横浜市における救急出場件数に対する関数動的予測モデリング横浜市立大学 医学部

三角 俊裕

横浜市立大学 医学部

三枝 祐輔

横浜市立大学 医学部

原 みゆひ

横浜市消防局

藤田 豊

横浜市消防局

吉塲 亮介

1. 概要
図１は，2008 年から 2019 年に記録された横浜市における救急車の出場要請データを示す．救急車の時間別出場件数は，午前中にピークを迎え，午後から夜間にかけて減少する特徴を持つ．また，休日・
休日明け・平日に分類すると，特に休日明けの午前中に増加する傾向にある．さらに，出場件数は夏と冬に比較的多く，春と秋に少ない季節変動性の特徴がある．2008 年以降は増加傾向にあり，今後もこの傾向が継続する場合，救急車の現場到着の遅延などの問題が危惧される．そのため，限られたリソースを効率よく運用する施策が求められている．

40

50

Weekday

30

30

40

50

Day after holiday

20
0

10

20
10
0

0

0

10

10

20

20

Count

Count

30

30

40

40

50

50

Holiday

0

5

2014

Date

2016

2018

5

10

20

0

5

10

20

700
600

650

2019
2018
2017
2016
2015
2014
2013
2012
2011
2010
2009
2008

550

35
30
25

450

15
2012

0

2019
2018
2017
2016
2015
2014
2013
2012
2011
2010
2009
2008

350

400

10
2010

20

Hour

5

300
2008

10

20

20

Mean counts

700
600
500
400

Total counts per day

15

500

10
Hour

Mean couns

5

800

0

5

10

Hour

15

20

2

4

6

8

10

12

Month

図 1: 左上：2008-2019 年からランダムに選択した 500 日分の時間別救急出場件数のスパゲッティプロッ
トと平均曲線，右上：2008-2019 年からランダムに選択した 500 日分の時間別出場件数を休日・休日明け・平日に分類して作成したスパゲッティプロットと平均曲線，左下：1 日あたりの合計救急出場件数の時系列プロット，中央下：時間別救急出場件数の年平均曲線，右下：月別救急出場件数の年平均曲線．

本研究では，一日のある時間までに観測された救急出場件数に基づいて，以降の出場件数予測を逐次的に行う動的予測法を提案する．動的予測法の文脈では，交通量予測 (Chiou, 2012) や子供の成長予測(Leroux, 2018) など，関数データ解析手法に基づくアプローチが注目されている．そこで，様々な特徴を持つ横浜市の救急出場データに対して，Leroux (2018) により提案された関数同時回帰モデルを拡張させた関数動的予測モデリング手法を提案した．提案手法は，横浜市の救急出場データに適用した．また，数値シミュレーションにより提案手法の有用性を検討した．

2. 関数動的予測モデル
n 日分の観測データを {ξij , xiq , zijr ; i = 1, . . . , n, j = 1, . . . , J, q = 1, . . . , Q, r = 1, . . . , R} とする．ここで，ξij = ξi (tij ) は第 i 日目の tij 時における出場件数，xiq および zijr = zir (tij ) は，それぞれ第 i 日目の q 番目の非時間依存共変量，tij 時における r 番目の時間依存共変量とする．また，出場件数の長期傾向と季節変動を考慮するために，第 i 日が属する年，月をそれぞれ tyik (k = 1, . . . , K)，tm
il (l = 1, . . . , L)とし，当該出場件数を ξikl (t) と定義する．このデータに対して，以下のモデルを仮定する．

ξikl (tij ) = f0 (tij ) +

Q
∑
q=1

xiq f1q (tij ) +R∑zir (tij )f2r (tij ) + g(tyik , tmil ) + bi (tij ) + εi (tij ).

(1)

r=1

ここで，f0 (·) は母平均関数，f1q (·) は q 番目の非時間依存共変量の係数関数，f2r (·) は r 番目の時間依
存共変量の係数関数，g(·, ·) は長期傾向と季節変動を表す関数，bi (·) は日毎の変動を表す変量効果関数，

εi (tij ) = εij は観測誤差である．また，bi (t) は共分散関数 C(s, t) = cov{(bi (s), bi (t)} をもつガウス過程，
εij は N (0, σ 2 ) にそれぞれ従い，bi (t) と εi (t) は互いに独立であると仮定する．さらに，時間依存共変量は，m
季節変動の特徴をもつ気温や湿度等を想定し，zirk (tij ) = γr (tij )+gr (tik)+bir (tij )+εir (tij ) のようなモデ

ルを (1) 式と同時に仮定する．ここで，γr (·) は母平均関数，gr (·) は季節変動を表す関数，bir (·) は日毎の変動を表す変量効果関数，εir (tij ) = εijr は観測誤差とする．bir (t) は共分散関数 Cr (s, t) = cov{bir (s), bir (t)}をもつガウス過程，εijr は N (0, σr2 ) にそれぞれ従い，bir (t) と εir (t) は互いに独立であると仮定する．
出場件数の長期傾向と季節変動は，g(ty , tm ) = ψtrend (ty ) + ψseasonal (tm ) + ψinteraction (ty , tm ) と表し，

ψtrend (·) は長期変動項，ψseasonal (·) は季節変動項，ψinteraction (·, ·) は交互作用項とする．時間別出場件
数の動的予測は，Leroux et al. (2018) を参考に最良線形不偏予測量を用いて行った．
モデリングの詳細，数値シミュレーション結果および実データへの適用結果は，当日報告した．
参考文献

Leroux, A., Xiao, L., Crainiceanu, C., and Checkley, W. (2018). Dynamic prediction in functional
concurrent regression with an application to child growth. Statistics in Medicine, 37(8):1376–1388.
Chiou, J.-M. (2012). Dynamical functional prediction and classiﬁcation, with application to traﬃc
ﬂow prediction. The Annals of Applied Statistics, 6(4):1588–1614.

概念の複雑性・信頼性予測のための受け入れ検査 (デジタル言語学)
得丸久文(カラハリプロジェクト)

1. デジタル言語学（生命進化の延長線上にある言語的知能）
1.1 原核生物からインターネットまで一貫する生命のデジタル進化(物理層➔論理層)
デジタルは 01(ON-OFF)の電気信号にもとづく 2 進計算に限らない．むしろ生命体が 4 種の核酸で紡いだ遺伝情報を自律分散的に進化させる生命進化こそがデジタルであり，2 進計算機と離散信号である音節を使った言語と知能の進化は生命進化の派生として位置づけられる．
表 1 は真核生物，脊椎動物，哺乳類，人類の進化を概観する．すべて進化は，よりよく生きよ
うとする生命体の意志が生みだす．第一段階は物理進化であり，小惑星衝突などによる生息環境の悪化が，それを乗り切るための低雑音環境(核，脳室，子宮)を生みだす．第二進化は論理進化
であり，環境が復元されると，低雑音環境(=高 S/N 比)での高ダイナミックレンジを利用した，
複雑で繊細な信号処理が可能となり，目に見えない論理進化が生まれる．

時期

誕生

物理進化

物理進化現象論理進化現象

2BA

真核生物

核，共生

核の保護と細胞質共生

転写後修飾，多細胞化

530MA

脊椎動物

脳室

視覚・聴覚

脊髄反射(脳室内免疫 NW)

66MA

哺乳類

子宮

保育器，授乳大きな脳，音声信号交信

130KA

現生人類

洞窟居住

ヒトの晩成化ヒトの真社会性化，鳴き声共有

66KA

言語的人類

音節

無限の語彙数，母音文法(母語の片耳聴覚)

5KA

文明人

文字

文明(時空を超えた共有)

概念(反射から群への進化)

50A

知恵ある人

Bit

ネットワーク検索

前方誤り訂正表 1 生命のデジタル進化概要(B=10 億，M=100 万，K=千，A=年前)
1.2 離散信号音節を活用するためには低雑音環境が必要であった
ヒトの言語は，真社会性動物が共同体内部で共有する特殊な鳴き声が，相互に周波数離散的な音素性を獲得し，母音アクセントによる時間的離散性をもつ音節を獲得して，無意識の文法処理が始まった．その後，消えない音節である文字と対話する音節である bit が生まれた後，その特徴を活かした論理進化を確立させ，知恵ある人の遺伝子に書き込む日がまもなく訪れる．
文法処理が無意識に自然に行われるのは言語学上の最大の謎であった．筆者の仮説は①
し手は，文法語の音節を識別子なしに挿入して修飾・接続を示す，②

話

聞き手は，母語を片耳で聴き，脳幹聴覚神経核の方向定位能力で文法語をベクトル解析する，というもの．(得丸 2014)
片耳聴覚への切り替えが生得的な遺伝子情報に基づいている(遺伝子に符号化されている)こと，熱帯雨林で生活するピダハンはその遺伝子が発現せず母語を両耳で聞くことから推論すると，文法のためには開墾地や都市など低雑音環境が必要である．概念は世俗から隔絶された僧院で生まれた．前方誤り訂正のためには，より低雑音の環境とより低速の反復処理が求められる．

2. 概念の複雑性を予測するために
2.1 言葉を処理する生物構造は脳室内免疫細胞ネットワークである
ヒトの言語的知能の基盤は，脊椎動物の脊髄反射回路であり，それは脳室内免疫細胞ネットワークである．五官の記憶はマイクログリア細胞に記銘されて大脳皮質に着床し，言葉の記憶とネットワーク記憶は脳脊髄液中を浮遊する B リンパ球(B 細胞)が保持する．

2.2 言葉の意味は 2 種類の記憶と 3 種類の論理によってフラクタルに複雑化する
意味の複雑化は，言葉と意味を二元的に結びつける B 細胞の内部論理が，1 対 1 の反射論理から，1 対全の群論理，全対全のネットワーク論理へとフラクタルに複雑化して実現する．言葉には 6(=2x3)通りの複雑性があることを聞き手は知っておく必要がある．
五官記憶(グリア細胞)

意味記憶(B 細胞，ネットワーク)

1 対 1(反射論理)

偶発的体験

現象の恣意的な命名

1 対全(群論理)

日常概念

科学的概念

全対全(ネットワーク論理)

日常概念操作

学際的な科学概念

3. 概念の信頼性は予測不可能．知識再生的に受容する
3.1 脊髄反射と真社会性の制約により，知能は誤り訂正しにくい
文字や bit 化された言語情報は，通信回線上で単なる文字列・bit 列となるため，改ざんが容易である．著者没後の偽書を見破ることは難しい．
真社会性動物であるため，師が弟子に知識伝達する場合，弟子は鵜呑みにするほかない．
脊髄反射回路で言語処理するため，一度誤った知識を獲得するとそれが基準になり，正しい情報を受け入れられなくなる．

3.2 知識再生的な受容のための受け入れ時のアルゴリズムを実装する
受け入れる知識は，話し手が自分で考えたことか(A)，他者から仕入れた知識か(B)A ➔ どのように考えたかの筋道を明らかにするB ➔ 誰が考えたのかを明らかにする➔ 自分でも同じように考えてみる➔ 知識・概念の発見者を見つける ➔A に戻る

この知識再生的受容アルゴリズムが遺伝子に書き込まれたとき，ホモサピエンスが誕生する．

参考文献：得丸(2014) 母語のモノラル聴覚と文法処理－例外としてのピダハン, 情処学会 2014-NL-219(23) PP 1-19
得丸(20H00576)概念の複雑性、ならびに『きわめて大量だが信頼性の保証がない言語情報』の取り扱いについて得丸(20H00576) 概念は群であるー学習における前方誤り訂正の必要性

予測のための最大エントロピーと最小ダイバージェンス
江口 真透
統計数理研究所
一般化線形モデル (GLM) は指数型分布族の中で標準パラメータから回帰関数へのリンク関数によって自然な線形モデルを導入している．更に混合効果モデル，一般化推定方程式などの工夫から広くデータ解析に用いられている．その理論的背景には最大エントロピーモデルと KL ダイバージェンスの関係の巧妙な利用がある．この発表では回帰の内容で，一般の U -エントロピーと U -ダイバージェンスを考えて最大エントロピーと最小ダイバージェンスの双対的な関係を追求する．これによって GLM の拡張を図り，GLM をより柔軟にする方法について考えたい．

1. はじめに特徴ベクトル X から目的変数 Y の予測問題を考える．X = x を与えたとき Y の条件付き密度 (質量) 関数の全体∫を P = {p(y|x)} と書くとき，回帰関数の空間 R = {µ(x) = yp(y|x)dΛ(y) : p ∈ P} が導かれる．GLM の枠組
みでは，分散パラメータを取り込まない形であれば

p(y|x) = exp{yθ(x) − ψ(θ(x))}と仮定される．cf. Hastie et al. (2009). 標準パラメータ集合が実軸 R まで拡張できると仮定すると標準リンク
モデル θ(x, α) = α⊤ x が無理のない形で回帰モデルが定式化できる．これより，標準パラメータから平均値パラメータの変換から回帰関数は µ(x, α) = η(α⊤ x) と表される．ここで η(θ) = (∂/∂θ)ψ(θ) とする．このように Yがベルヌイ分布，ポアソン分布，ガンマ分布などに従っているときでも，θ(x) のモデル化によって自然な回帰モデルが提案されている．GLM の枠組みにおいて (1) はエントロピー最大分布となっており，パラメータ推定は尤
度法，あるいはクルバック・ライブラー（KL）ダイバージェンス最小化が実行されている．この関係を一般化して，そのクラスの中でデータをより柔軟に学習する方法を求めたい．ここで P 上のエントロピーは∫
H0 (p) = p(y|x) log p(y|x)dydx
で与えれ，対数尤度関数はデータ D = {(xi , yi )ni=1 } に対してℓ(α, D) =n∑{yi α⊤ xi − ψ(θ(α⊤ xi ))}i=1と書けて yi に関して線形であることが特徴となっている．

2. 一般化エントロピーとダイバージェンス
エントロピーとダイバージェンスの一般化のため生成関数の集合を U = {U : R → [0, ∞) : U ′ (t) > 0, U ′′ (t) > 0}
用意する, cf. Eguchi (2006). スカラー倍の不定性を除くために U ′ (1) = 1 を課す．U (t) を単調増加な凸関数とする．この U の U を使って，空間 P に交差エントロピーを∫∫
HU (p, q) = − {p(y|x)ξ(q(y|x)) − U (ξ(q(y|x))}p(x)dΛ(x, y)と定めて，U -交差エントロピーと呼ぶ．ここで ξ(s) は U ′ (t) の逆関数とする．このとき，

HU (p, q) ≥ HU (p)が成立し，等号は p = q のときに限ることが知られている．ここで HU (p) = HU (p, p) とし，HU (p) を U -エント
ロピーと呼ぶ．このことから U -ダイバージェンス

DU (p, q) = HU (p, q) − HU (p)が定まる．与えられたデータ D = {(xi , yi )ni=1 } から統計モデル M = {p(y|x, θ) : θ ∈ Θ} が建てられたとき，U 損失関数を∫n∑
LU (θ, D) = −{ξ(p(yi |xi , θ)) − U (ξ(p(y|xi , θ)))dΛ(y)}i=1と定めて未知パラメータ θ の U -推定量 θ̂ U = argmaxθ∈Θ LU (θ, D) を導入する．U -損失関数 LU (θ, D) は交差エントロピー HU (p, q) の p(y|x) はデータ D の経験分布関数，q(y, x) はモデル p(y|x, θ) を代入したものとなっている．
U の中から U0 = exp を取れば HU0 (p) はボルツマン・シャノンのエントロピー，DU0 (p, q) は KL ダイバージェ
ンス，LU (θ, D) は負の対数尤度関数に帰着され，GLM の基礎をあたえる．一方で，U の中から
β+1

1
(1 + βt) βUβ (t) = β+1を採用すると，ξ(s) = β1 (sβ − 1) となるので，べきクロスエントロピーは
∫
{1
β
β+1
1
Hγ (p, q) = −}p(x)dΛ(x, y)β p(y|x)[q(y|x) − 1] − β+1 q(y|x)
が導かれる．
このように生成関数集合 U に対してモデルと推定の組を U × U から (U0 , U1 ) を取り，U0 -エントロピー最大モ
デル MU0 と U1 -推定量 α̂U1 を取ると対角集合 ∆ := {(U, U ) : U ∈ U } に共通な構造が現れる．すなわち，U0 エントロピー最大モデルは()
MU0 = {p0 (y|x, α) := U0′ yα⊤ x − ψ0 (α⊤ x) : α ∈ Rd }で与えられ，U1 -ダイバージェンス推定の損失関数はモデル M = {p(y|x, α) : α ∈ Rd } に対して∫n∑
LU1 (α, D) = −{U1′ −1 (p(yi |xi , α)) − U1 (p(y|xi , α))dy}i=1で与えられる．従って，(U0 , U1 ) = (U, U ) と取れば (1) はLU (α, D) = −
⊤

⊤

∫

n
∑

{yi α⊤ xi − κ(α⊤ xi )}i=1と書ける．ここで κ(α x) = ψ0 (α x) − U (yα⊤ x − ψ0 (α⊤ x))dy. このように LU (α, D) は負の対数尤度関数と等価な関係が成立する．対角集合 ∆ の中で最適なものを探す方法についてはいくつか考えられるが決定的なものは未だ得られてない．応用として生物種の生息分布の推定のためのマックスエントを β-マックスエントに拡張されたものが考えられている, cf. Komori & Eguchi (2019)．今後，様々な応用が成される中で有効なものが提案されることが期待される．

参考文献
Eguchi, S. (2006). Information geometry and statistical pattern recognition. Sugaku
Expositions, 19(2), 197-216.
Komori, O., & Eguchi, S. (2019). Statistical methods for imbalanced data in ecologicaland biological studies. Springer Japan.

臨床研究における統計家の役割
静岡県⽴静岡がんセンター
臨床研究⽀援センター
野津昭⽂
臨床研究にはさまざまな種類(表 1)がある。診療記録を⾒返して、治療の効果や安全性について集計を⾏う研究は、後ろ向きの研究と⾔われ、臨床研究では重要な役割を担っている。
しかし、後ろ向き研究には、⼀般には多くのバイアスが含まれており、結果の信頼性は⾼いとは⾔えない。そこで、収集すべき項⽬や評価⽅法をあらかじめ決めてから、データを収集する研究があり、前向きの研究と⾔われる。後ろ向き研究とは異なり、前向き研究では、結果として得られたデータの質を担保するために、多くの⼯夫がなされている。例えば、医療従事者だけでなく、データの管理を⾏うデータマネージャー、研究の進捗を管理する研究事
務局、データ収集の補助を⾏う臨床研究コーディネーター、統計解析を⾏う統計家が関わることがある。これらの職種がそれぞれの領域で専⾨的な知識を持って関わることで、臨床研究によって質の⾼いデータを得ることができる。統計家は、研究の計画、実⾏、解析、解釈のすべてに少なからず関わることになる。
予測モデルの作成についても、上述のように、どのようにデータを取得するかという問題がある。後ろ向きの研究から得られたデータを使う場合もあれば、前向きの研究から得られたデータを使う場合もある。予測モデルの作成は、ともすれば、統計モデルの話に限定して語られることがあるが、どのように取得されたデータを、どのように使⽤するかという観点から、すなわち研究デザインの観点から議論することも重要である。統計家は、研究⽬的を果たすために、必要であれば研究組織についても考える場合もある。つまり、ここでも、研究の計画、実⾏、解析、解釈のすべてに少なからず関わることになる。
本発表では、病院で勤務している統計家が、臨床研究の試験デザインや臨床研究を進める際の組織、さらにそこでの統計家の役割について、当院の例を通して話をした。臨床研究は段階を踏んで進められること、各段階ですべきこと・⾏わないことについても解説した。さ
らに、病院の医療従事者と予測モデルを作成する研究を⾏った際に感じたことや、当院以外の事例も参照しつつ、研究デザインと予測モデルの関連についても話をした。
表 1：臨床研究の種類

観察研究
後向き研究

介⼊研究

症例報告、症例集積報告、
コホート研究、症例対照研究

前向き研究 コホート研究、症例対照研究

パイロット研究、忍容性試験、第I, II, III, IV相試験参考⽂献
井村裕夫(監修): NIH 臨床研究の基本と実際, 丸善, 2016
⽊原雅⼦/⽊原正博(訳): 医学的研究のデザイン, メディカル・サイエンス・インターナショ
ナル, 2014
⼤橋靖雄 荒川義弘: 臨床試験の進め⽅, 南江堂, 2006

ガンマクレームによる資産過程における破産確率の漸近近似
千葉大・融合理工学府 岡山 大成
千葉大・理学研究院

内藤 貫太

資産過程のクラメル・ルンドバーグモデルでは，クレーム額を IID 確率変数とし，時刻 t までのクレーム
回数を強度 λ のポアソン過程 N = (Nt )t≥0 を用いて表現する．ここで，i 番目のクレーム額を表す確率変数

Yi の分布を FY とし，E[Yi ] = µ とする．初期資産を x ≥ 0，安全付加率を θ とすると，保険会社の資産過程
X = (Xt )t≥0 は
Xt = x + (1 + θ)λµt −

Nt
∑

Yi

i=1

と表される．初期資産 x の Xt についての破産時刻を τ (x) = inf{t > 0 : Xt < 0} とし，有限時間破産確率を

ψ(x, T ) = P(τ (x) ≤ T ) とする．このモデルが純益条件を満たすとき，破産確率 ψ(x) = limT →∞ ψ(x, T ) のラ
プラス変換はLψ(s) =µθ
1−
s sµ(1 + θ) − 1 + mY (−s)

となることが知られている (清水 (2018))．ここで mY (s) は Y のモーメント母関数である．
資産過程 Xt において，自然数 n に対して，

Yiθθ → √ , λ → nλ, Yi → √nnという置き換えによって得られる資産過程

(
)
Ñt
θ
1 ∑
µ
X̃t = x + 1 + √Yi(nλ) √ t − √
n
n
n i=1
をスケーリングによる資産過程と呼ぶ (Cohen and Young(2020))．ここで，Ñt は強度 nλ のポアソン過程であり，n は例えば観測期間に対応している．極限定理として，dX̃t −→ x + Wt (n → ∞)が成り立つことが知られている．ここで，Wt は (θλµ, λE[Y 2 ])-ブラウン運動である (Schmidli(2017,Proposition

5.19), 清水 (2018,6.3.2 節))．
n に依存している X̃t の破産確率を ψn (x) とする．n = 1 のとき，X̃t は元の資産過程 Xt そのものである．したがって，Xt の破産確率 ψ(x) は，n = 1 のときの X̃t の破産確率 ψ1 (x) である．いま，極限定理が成り立って
おり，この極限資産過程 (ブラウン運動) の破産確率を ψD (x) とすると，形式的に ψn (x) → ψD (x) (n → ∞)である．よって，ψ(x) = ψ1 (x) が ψD (x) で近似できているのかが問題となる．ψn (x) → ψD (x) (n → ∞) であり，ψ(x)=ψ1 (x) であるから，ψ(x) の ψD (x) による近似は ψn (x) を陽に求めることで考察できる．

本発表では特に，形状パラメータが α，尺度パラメータが β のガンマ分布にクレーム額が従う場合を考える．
そのとき，n → ∞ において ψn (x) は以下の漸近展開をもつ：

{

}
( )
1
1
1
ψn (x) = ψD (x) 1− √ A1 (x|θ, α, β)+ A2 (x|θ, α, β) +o
n
n
n
ここで，

2θβ
ψD (x) = e−γx , γ =
,
α+1
{
}
2θ(α + 2) θ(α − 1)
A1 (x|θ, α, β) =+1(x = 0) (1 − γx),3(α + 1)3(α + 1){}θ2 (α + 2)3(α2 − 1)
2 2
A2 (x|θ, α, β) =
2(α
+
2)γ
x
−
3(3α
+
5)γx
+
6(α
+
1)
+
1(x
=
0)
9(α + 1)2
α+2である．この展開は Cohen and Young(2020) で議論された展開の一つの一般化となっている．
また，ψn (x) について，任意の x ≥ 0 において

ℓn (x) < ψn (x) < un (x)が成り立つような ℓn (x) と un (x) を与えることができる (Cohen and Young(2020))．ψn (x) の漸近展開式において，o(n−1/2 ) 項を落とした{}1
ψn,1 (x) = ψD (x) 1 − √ A1 (x|θ, α, β)nに対して，任意の x ≥ 0 において

ℓn (x) < ψn,1 (x) < un (x)が成り立つ．これによって，|ψn (x) − ψn,1 (x)| < 1(ψn (x) ≥ ψn,1 (x))|ψn (x) − ℓn (x)| + 1(ψn (x) < ψn,1 (x))|ψn (x) − un (x)|が成り立ち，ψn,1 (x) は Cohen and Young(2020) で与えられた限界よりも ψn (x) の良い近似を与えている．
本講演において，Ai (x|θ, α, β)(i = 1, 2) の導出について解説を与えるとともに，この展開式の精度と応用例
について報告を行った．

参考文献[1] 清水泰隆 (2018). 保険数理と統計的方法，共立出版．
[2] Cohen, A. and Young, V. R. (2020). Rate of convergence of the probability of ruin in the CramérLundberg model to its diﬀusion approximation, Insurance: Mathematics and Economics, 93, 333-340．
[3] Schmidli, H. (2017). Risk Theory. Springer Actuarial Lecture Notes, Springer.

⽣物分布の時空間予測：機構論的モデル構築に向けた課題
楠本聞太郎（九州⼤学農学研究院）
将来的な⽣物多様性損失を防ぐための保全戦略の構築と実践が急務となっている。⽣物多様性は地理的空間上に不均等に分布しており、我々は限られた保全リソースの最適配分を考えなければならない。⽣物の空間分布情報は、⽣物多様性のホットスポットを特定し、様々な社会経済的制約の下で、最⼤限の保全効果を得るための空間プランニング
にとって必須の情報である。保全効果を⻑期的に持続させるためには、現在の⽣物の分布だけでなく、気候変動による⽣息適地の変化も⾒越した保全エフォート配分が重要になる。これを実現するためには、それぞれの種が気候変動にどのように応答し、その結果として⽣物多様性の分布がどのように変化するかを予測できなければならない。本講
演では、種分布モデリングによる時空間予測のトレンドをレビューし、現状の情報ギャップを指摘する。そして、情報ギャップを緩和するための解析的な⽅法と、データ収集の取り組みを紹介した。
現在、最も広く普及している分布モデリング⼿法は、⽣物の観測（在）データとそれが得られた場所の環境要因との相関関係に基づくモデリングである。相関ベースの⽅法は、特別な⽣物学的メカニズムを考慮しないにも関わらず、多くの場合で、対象⽣物の現在の空間分布（在・不在）をよく判別する。この⽅法は、最⼩限の情報（在データと環境データ）でモデリングできることから、多数の⽣物種を扱う研究にも適⽤しやすいのが強みである。⼀⽅で、⽣物学的メカニズムを考慮しないことへの批判もある。相関ベー
スモデルに対して、機構論的モデルは、⽣物の分布レンジを規定するメカニズムとそれらを制御する環境条件の関係を明らかにし、地理的空間に投影する。実際の地理分布は、複数のメカニズムが複合的に働いた結果であり、全てのメカニズムを統合したモデルを構築することは困難である。よって、純粋な機構論的モデルは、部分的にしか分布の制限要因を反映せず、分布レンジを過⼤評価する傾向がある。このように、相関ベースモ
デル、機構論的モデルいずれにも制約や問題点がある。2 つのアプローチの予測結果の違いは、構築したモデルを別の地理的空間や、将来環境に投影するときに顕著になるようである。このようなモデル間の不⼀致性について、どちらのモデルが正しいかを判断するのは難しい。機構論的モデルには⽣物学的なリアリズムがあるが、メカニズムに対する仮定が間違っていたり、重要な⽣物学的プロセスを⾒落としていたりすると、相関ベースよりもむしろ悪い予測結果になることもある。モデル間の不⼀致性は、モデリン
グに⽤いた環境要因では捕捉できないメカニズムが働いていることを⽰唆しているため、異なるアプローチで構築されたモデルを⽐較しながら、⾒落とされている⽣物学的メカニズムを考察する必要がある。
双⽅の強みを活かすために、相関ベースモデルに⽣物学的なリアリズムを反映させたハ
イブリッド型アプローチも提案されている。ハイブリッド型アプローチの最も単純な例
は、メカニズムを直接的に制御する環境変数を選択する⽅法である。そのほか、分布を制限するメカニズムの情報を使って、相関ベースモデルの学習を⾏う領域を操作する⽅
法や、局所スケールで得られた個体群レベルの⾼解像度データを使って相関ベースモデ
ルを構築し、それを広い時空間スケールに投影する⽅法が提案されている。
⽣物学的メカニズムに関するデータを得るには、多くの実験・観測情報が必要であり、ほとんどの分類群で、そのような情報は利⽤可能でない。このことが、機構論的モデル
やハイブリッドモデルの多種系への適⽤への障壁になっている。これに対する解析的な解決策として、階層モデルやインバースモデルによって、異なるタイプのデータをフュ
ージョンする⽅法が着⽬されている。階層モデルの例としては、多種系で利⽤可能な種
特性や系統情報を⽤いて、メカニズムに関連するニッチパラメータを補完（推定）する
⽅法や、環境応答の種間差を潜在変数として多種の分布予測を⾏う連結種分布モデル（joint-SDM）がある。インバースモデルは、「⽣物学的メカニズムに関するデータ→機
構論的モデル→地理的空間への投影」という流れとは逆に、「地理分布→機構論的モデ
ルのパラメータ推定」という流れで解析を⾏う。これは、少数の種についての詳細な観測・実験データと、多種系で利⽤可能な在データの両⽅の強みを活かせる点で有望である。
最後に、情報不⾜に対する究極の解決策は、新しいデータの取得である。本講演では、
⽇本産樹⽊種の種⼦に関する⽣理学的情報を収集し、それに基づいて機構論的（ハイブ
リッド型）分布モデルの構築を⽬指すプロジェクトを紹介した。樹⽊のような固着⽣物の場合、ある場所での⽣存個体の在情報が分布モデリングに⽤いられる。しかし、その場所の環境が、種⼦の貯蔵・発芽から、個体の成⻑・⽣存、繁殖フェノロジーまで、統合的に適した環境なのか、将来的にも好適環境が保たれるかどうかはわからない。本研
究では、既存の⽂献情報、データベースに加え、圃場・実験室での発芽実験を通じた、種⼦の温度依存的な特性に関する情報収集と、それらと分布情報との関係の検証を⾏う
予定である。

Some recent topics for Preferential Sampling
Shinichiro Shirota
October 23, 2021
The basic issue of preferential sampling (PS Diggle et al., 2010) is bias in the samplingof spatial locations where point referenced response data are collected and the potential
impact on inference regarding the response surface over the region. The canonical illustrative
example addresses the objective of inferring about environmental exposures. If environmental
monitors are only placed in locations where environmental levels tend to be high, theninterpolation based upon observations from these stations will necessarily produce only highpredictions. The obvious remedy lies in suitable spatial design of the locations. For example,
a random or space-filling design (Saltzman and Nychka, 1998) for locations over the regionof interest is expected to preclude such bias.
However, sampling may not be designed in this fashion. Environmental researchers may
install monitoring stations where they expect to find high exposure levels; ecologists may
tend to sample where they expect to find individuals. Recognizing the possibility of such bias,
can prediction be revised to adjust for it? This is the intention of PS modeling. Specifically,
there are typically two questions: (i) is there evidence of a PS effect? and (ii) can we improveprediction in the presence of PS? That is, we do not seek to remove PS, we do not propose
to revise the data collection. Rather, we seek to acknowledge its presence and attempt tomitigate its impact.
We introduce the Bayesian modeling approach for PS based on the point processes (especially, log Gaussian Cox processes, Møller et al., 1998) for observed locations, then briefly
review some recent topics on PS. Our main contribution here is to expand the issue of sampling bias to the case of bivariate response. If responses at location s, say (Y1 (s), Y2 (s)) aredependent then prediction of say Y1 (s0 ) at an unobserved location s0 can benefit from using
all of the observed Y2 ’s in addition to all of the observed Y1 ’s.
Then, the question here becomes the following. Suppose Y1 and Y2 are strongly dependent, that is, for a given region D, realizations of the response surfaces, Y1 (s) : s ∈ D andY2 (s) : s ∈ D are highly correlated. Then, if there is sampling bias in the subsets, S1 andS2 in D where Y1 and Y2 are sampled, respectively, can we demonstrate the presence of a1bivariate preferential sampling effect? Can we improve co-kriging through a model which
captures bivariate preferential sampling?
We consider two illustrative examples here. One focuses on realizations over a spatial region of ozone and temperature in the summer season (May-September). It is well establishedthat high temperatures are a driver of ozone formation, suggesting expected strong positive
dependence. The other considers, for a site, mean diameter at breast height (MDBH) overthe site and trees per hectare (TPH) associated with the site. Here, the “constant yield law”
(Weiner and Freckleton, 2010) argues that the total yield/biomass on a plot at equilibrium, isroughly constant regardless of the number of individuals, i.e., size of individuals will decrease
as density increases, suggesting strong negative dependence.
In the bivariate setting we have Y1 with associated S1 and Y2 with associated S2 . This
raises two possibilities: (i) S1 = S2 and (ii) S1 6= S2 . Above, the ozone and temperature
fall under (i) since they are collected at the same set of monitoring stations. We imagine asingle point pattern but bivariate geostatistical response. The MDBH and TPH fall under(ii) because the sites for the measurements can be different due to different data collection
protocols.

References
Diggle, P., R. Menezes, and T. Su (2010). Geostatistical inference under preferential sampling. Journal of the Royal Statistical Society, Series C 59, 191–232.
Møller, J., A. R. Syversveen, and R. P. Waagepetersen (1998). Log Gaussian Cox processes.
Scandinavian Journal of Statistics 25, 451–482.
Saltzman, N. and D. Nychka (1998). DI, a design interface for constructing and analyzing
spatial designs. Springer Science and Media.
Weiner, J. and R. P. Freckleton (2010). Constatn final yield. Annual Reviews of Ecology,
Evolution, and Systematics 41, 173–192.

2

近似ベイズ計算による野生動物の集団サイズ推定
生成過程が複雑なデータにどう対処するか？
矢島 豪太・中島啓裕（日本大学生物資源科学部）
はじめに
野生動物の多くの種は，複数個体から成る集団を形成して生活している．その集団サイズ（単位集団内の個体数）は，動物の社会生態を明らかにするうえで基礎となる情報である．
また，個体数密度を推定したり，個体群の動態を予測したりするうえでも，集団サイズとその時間変化についての情報は不可欠である．しかし，観測時に集団内のすべての個体を数え上げることは困難であるにもかかわらず，従来の研究の多くは，こうした「偽陰性」を考慮
せず，観測値をそのまま真の集団サイズと見なして解析・推論を行ってきた．
近年，生態学の分野で自動撮影カメラが広く利用されるようになってきた（Burton et al.
2015）．自動撮影カメラは，赤外線センサーに動物が反応すると自動で写真を撮影するシステ
ムのことである．自動撮影カメラは一度設置してしまえば，昼夜問わず記録を行えること，
また動物に警戒されずに記録を行えることから，野生動物のデータ取得における大きなブレークスルーとなった．
しかし，自動撮影カメラによる観測データから集団サイズを推定することは，従来の直接観察よりもいっそう困難である．第一に，撮影可能な面積が限られているため，個体の検出
率が群れの広がりの影響を受け，個体の検出・非検出が独立に起こらない．第二に，陸生の社会性動物の多くは，互いに重複のある固定的な行動圏を持つことが多く，自動撮影カメラのような受動型の検出機器には，複数集団が複数回撮影されるのが普通である．このため，
集団サイズの推定には，その地点に生息する集団数と集団サイズを同時に推定する必要が生
じる．すなわち，カメラのデータを利用して，正確な集団サイズの推定を行うためには，集団内の個体が観測されるプロセスを明示的に表現した新たな統計モデルが必要である．さらに，複雑なモデルのパラメータ推定を効率よく行うためのアルゴリズムの検討が不可欠である．
そこで，本研究の目的は，(1) 自動撮影カメラで取得したデータに適用可能な，集団サイ
ズ推定のための統計モデルを構築すること，(2) 考案したモデルを房総半島（千葉県南部）のイノシシ個体群に実際に適用することの 2 つである．モデルのパラメータ推定は，近似ベイズ計算（Approximate Bayesian Computation，以下 ABC）を利用した．
統計モデルと推定アルゴリズム
統計モデルの構築は，真の集団数とサイズが決定される生態学的なプロセス（生態プロセ
ス）と，
観測値を得るまでの観察プロセスを明示的に分離した階層ベイズモデルを利用した．
まず初めに生態プロセスを考える．調査エリア内に，空間的に独立な𝐼個の自動撮影カメラが設置されているとする．あるカメラ𝑖 𝑖 1, … , 𝐼 によって撮影される可能性のある（カメラと行動圏が重複する）集団数𝐾 は，パラメータ𝜇をもつポアソン分布に従うとする．また，各集団の真の集団サイズ𝑁 , 𝑘 1, … , 𝐾 も，パラメータ𝜆をもつポアソン分布に従うとす
る．
𝐾 ~ 𝑃𝑜𝑖𝑠𝑠𝑜𝑛 𝜇𝑁 , ~ 𝑃𝑜𝑖𝑠𝑠𝑜𝑛 𝜆
次に，観測プロセスを考える．𝑖番目のカメラの前を利用する𝑘番目の集団がカメラに観測される回数𝐶 , もやはりパラメータ𝜐をもつポアソン分布に従うと仮定する．従って，𝑖番目の∑カメラでは，𝐿𝐶 , 回の観測が得られることになる．また，観測集団サイズ𝑦 , 𝑗
1, … , 𝐿 は，個体の検出が独立ではないことを想定して，真の集団サイズ𝑁 , と𝛼，𝛽をパラメータにもち，二項指標変数𝑧 , を導入したベータ二項分布と見なす．
𝑦 , ~ 𝐵𝑖𝑛𝑜𝑚𝑖𝑎𝑙 𝑁 , , 𝑝 ,
𝑝,

𝑝′ ,

𝑧,

1.0

1

𝑧,

𝑝′ , ~ 𝐵𝑒𝑡𝑎 𝛼, 𝛽𝑧 , ~ 𝐵𝑒𝑟𝑛𝑜𝑢𝑙𝑙𝑖 Ω𝑧 , は，Ωをパラメータとするベルヌーイ分布に従うとする．この変数を導入したのは，野生1）と高い状態（𝑧 ,0）に二分できることが多い動物の集団は，集合性が低い状態（𝑧 ,と考えられたためである（例えば，採食時は個体間距離が離れやすいのに対し，移動時は１列になり集合性が高くなる）．ここでは，𝑧 ,0のときに検出率は 1 とした．すなわち集合
性が十分に高い場合は，カメラに全ての個体が撮影されると仮定した．
このモデルを既存の確率的プログラミング言語（Jags や Stan などの MCMC を利用するためのソフトウェア）で実装するためには工夫が必要になる．なぜなら，MCMC のステップ毎に集団数𝐾 の値が変化し，それに伴い真の集団サイズ𝑁 , の長さも変化するためである．そこで，より効率的な推定手法として，近似ベイズ計算（Approximate Bayesian Computation，
以下 ABC）を利用した．ABC は，適切な要約統計量を選択し，データ生成可能なモデルであれば，尤度関数を評価せずにパラメータの事後分布を推定できる（小林 2011）．ここでは推定効率を考え，逐次モンテカルロ法を応用した ABC–SMC（Toni et al. 2009）を利用することにし，観測集団サイズ𝑦 , と観測回数𝐿 の平均と分散を要約統計量として使用した．
当日の発表では，モンテカルロ・シミュレーションによるモデルの性能を検証した結果を報告した．また，千葉県房総半島で実際に撮影されたイノシシの集団に適用した結果についても報告した．
参考文献[1]. 小林弦矢．
「ABC サンプラーの最近の発展について」．日本統計学会誌，41，pp. 155–180，
2011．
[2]. Burton, A. C., Neilson, E., Moreira, D., Ladle, A., Steenweg, R., Fisher, J. T., Bayne, E. andBoutin, S. 2015. Wildlife camera trapping: a review and recommendations for linking surveys toecological processes. Journal of Applied Ecology 52: 675–685.
[3]. Toni, T., Welch, D., Strelkowa, N., Ipsen, A., and Stumpf, M. P. H. 2009. Approximate Bayesiancomputation scheme for parameter inference and model selection in dynamical systems. J. R. Soc.
Interface, 6, 187–202.

環境 DNA による生物多様性評価のための階層モデルとベイズ決定分析
深谷 肇一, 今藤 夏子, 松崎 慎一郎, 角谷 拓
国立環境研究所
fukaya.keiichi@nies.go.jp

環境 DNA（eDNA）メタバーコーディングは，非侵襲的かつ費用対効果の高い生物多様性の計測方法として広く応用されている．しかし，その多段階のワークフローの中では様々な要因で偽陰性が生じるため，

eDNA メタバーコーディングによる種の検出は完全ではない．eDNA メタバーコーディングの多段階のワークフローにおける不完全な検出はまた，研究デザインの問題とも関連する．具体的には，利用可能な資源を各段階にどう配分すれば調査効率を最適化できるかという問題である．
本研究では，対象地域内の複数のサイトで標本を収集する eDNA メタバーコーディング研究を対象に，多種サイト占有モデルの拡張を提案する (Fukaya et al. 2021)．二値の検出データを扱う従来の多種サ
イト占有モデル (Dorazio & Royle 2005) とは異なり，このモデルはハイスループットシーケンサーの出力である配列リード数の変動を記述するものである．提案モデルは eDNA メタバーコーディングの階層的なワークフローと種間の異質性を明示的に説明する部分モデルから構成された階層モデルであり（Fig. 1），
ワークフローの異なる段階における種の検出可能性の変動要因を分析可能である．また，ベイズ決定分析の枠組みを導入することで，限られた予算の下で種検出の有効性を最適化する研究デザインをモデルとそのパラメータの事後分布に基づいて特定できる．
具体的に，関心のある調査領域から J 個のサイトが選択され，eDNA メタバーコーディングを用いて合計 I 種の出現が調査される状況を想定する．サイト j (= 1, . . . , J) で得られた反復標本

k (= 1, . . . , Kj ) における各種の配列リード計数を yjk = (y1jk , . . . , yIjk )，シーケンス深度を∑
Njk = i yijk と表す．モデルは以下のように表される．
yjk ∼ Multinomial(Njk , πjk )uijk rijk
πjk = ∑I
m=1 umjk rmjk
rijk ∼ Gamma(ϕi , 1)(1)uijk ∼ Bernoulli(zij θi )(4)zijk ∼ Bernoulli(ψi )(5)(log ϕi , logit θi , logit ψi ) ∼ N (µ, Σ)

(6)

(2)
(3)

uijk と rijk はそれぞれ，サイト j の反復 k における種 i の DNA 配列の有無と，（配列が反復に含まれる場合
は）その相対頻度を表す潜在変数であり，zij はサイト j における種 i のサイト占有を表す潜在変数である．

1

Fig. 1. モデルの概要図．潜在変数 z（種のサイト占有を表す），潜在変数 u（反復標本における種の DNA
配列の有無を表す），データ y （配列リード数）の間に，種レベルのパラメータ（サイト占有確率 ψ ，配列
捕捉確率 θ，配列相対優占度 ϕ）で特徴付けられる階層的な依存構造を仮定する．配列リード数は，反復ごとの各種の配列の有無（u）と優占度（r）を条件とした多項分布に従うと仮定する．種レベルのパラメー
タは，平均ベクトル µ と共分散行列 Σ を持つ群集レベルの正規事前分布に従うと仮定する．Fukaya et al.
(2021)，CC ライセンス 表示 4.0 国際．
ϕi ，θi ，ψi はそれぞれ，種の配列の相対優占度，捕捉確率，および種のサイト占有確率と解釈でき，いずれもワークフローの異なる段階における検出可能性の変動を特徴付ける種固有のパラメータである（Fig. 1）．
このモデルを霞ヶ浦流域の淡水魚群集に適用したところ，種の検出可能性には顕著な不均一性があり，特定の種が偏って検出される恐れがあることが明らかとなった．サイト占有確率が低い種は捕捉確率が低く，配列の優占度も低いため，検出されにくい傾向にあった．期待される配列リード数は種によって最大 23.5 倍も異なると予測された．この結果に基づいて異なる研究デザインの間で種検出の有効性（期待される検出種数）を予測した結果，各反復に数万リードが確保されている限り，サイト内反復を複数確保することが種検出の有効性を高める上で好ましいことが示唆された．
本枠組みによって，eDNA メタバーコーディングがより誤差に対して頑健なものとなり，生物群集をより効率的にモニタリングできるようになることが期待される．

References
Dorazio, R.M. & Royle, J.A. (2005) Estimating size and composition of biological communities bymodeling the occurrence of species. Journal of the American Statistical Association, 100, 389–398.
Fukaya, K., Kondo, N.I., Matsuzaki, S.S. & Kadoya, T. (2021) Multispecies site occupancy modelingand study design for spatially replicated environmental DNA metabarcoding. Methods in Ecology and
Evolution. (to appear)2

Quantifying biodiversity change and its implication
Hideyasu SHIMADZU
Department of Mathematical Sciences, Loughborough University, UK
H.Shimadzu@lboro.ac.uk

Introduction
Quantifying biodiversity and its change has been a central interest in ecological sciences, and alarge number of metrics have been developed. The modern foundation of measuring biodiversity commonly consists of three ways:
α-diversity: the number of species found within a location;
β-diversity: the number of species commonly (not) found between locations; andγ-diversity: the number of species in a whole region.
Ideally, biodiversity metrics representing these three aspects should recognise the dimensionality of change and the ecological processes that lead to shifts in community composition.
However, it is not always the case and, in fact, even widely used measures such as the numberof species can fail to capture some changes (Shimadzu, 2018). This fact underlines an urgent
need for a better understanding of the nature of those metrics in terms of depicting biodiversity
changes.

Temporal changes in ecological communities
Let λ i ( t ) > 0 be the expected abundance of the i-th species at time t. The expected totalabundance of the community is then the sum of each species’ expected abundance as λ ( t ) Psi1 λ i ( t ) . Consider now an equation for the expected abundance of the i-th species as
λ i ( t + dt )  (1 + α i ( t ) dt ) λ i ( t ) ,(1)

where α i ( t ) is the instantaneous change rate at time t that drives the increase or the decrease of
species abundance. The change rate is a legitimate indicator to quantify, literally, the turnoverof the i-th species, based on its abundance, for the time increment dt. However, as the change
rate is in most cases unknown, it needs to be determined by the trajectory of the expected
abundance, λ i ( t ) > 0, itself as
α i ( t ) dt 

λ i ( t + dt ) − λ i ( t ) dλ i ( t ) d log ( λ i ( t )) .
λi (t )λi (t )

Although the continuous form of the expected abundances, λ i ( t ) , i  1, 2, . . . , s, is unknown,
we can estimate it using observations at discrete time points. We therefore rewrite Equation(1) in a discrete form representing a community dynamics model over a relatively short timeinterval, between times t and t + h (h > 0), within which the difference equation is still valid.
By integrating the instantaneous change rate, α i ( t ) , over the time interval, the community
dynamics can be described as
s
X

*1 +
λ (t + h ) 
i1 ,

Z t+h
t

!!
s
X
λi (t + h )+α i ( v ) dv λ i ( t ) 
1 + logλi (t ).
λi (t )i11Note that λ ( t )  si1 λ i ( t ) . The change rate and the instantaneous change rate may therefore
hold the following relationship for the short time period h,Psλ (t + h ) − λ (t ) X
λi (t + h )log
pi (t ) ,λ (t )λi (t )

!

(2)

i1

where p i ( t )  λ i ( t ) /λ ( t ) is the relative abundance of the i-th species at time t. The collectionof the relative abundances of the community — the relative abundance distribution p ( t ) P{p1 ( t ) , p 2 ( t ) , . . . , p s ( t ) } — satisfies the conditions that p i ( t ) > 0 and si1 p i ( t )  1. From
Equation (2), Shimadzu et al. (2015) define the turnover measure, D, between times t andu, ( u > t ) asD (t : u ) sXdi ( t : u ) 

i1

s
X


log

i1

 −

s
X
i1

λi (u )pi (t )λi (t )



(3)



pi (t )λ (u )log
p i ( t ) + log
pi (u )λ (t )!

 D1 (p ( t ) : p ( u )) + D2 ( λ ( t ) : λ ( u )) .
If the expected abundance, λ i ( t ) , is modelled via GLMs with the log link function aslog ( λ i ( t )) m
X

βi j x j (t ) ,j1

where {x j ( t ) } are environmental variables, and {β i j } are the parameters to be estimated, Equation (3) reveals that it becomes possible to identify drivers that influence the turnover measure,D, in an additive form as
D (t : u ) 

m X
s
X

β i j ( x j ( u ) − x j ( t )) p i ( t ) .

j1 i1

Since the contribution of the species and the environmental factors are all additive, puttingdi ( t : u ) mX

β i j ( x j ( u ) − x j ( t )) p i ( t ) andj1d j (t : u ) sX

β i j ( x j ( u ) − x j ( t )) p i ( t ) ,i1

the contribution of the i-species and of the j-th environment variable to the turnover measure,D, is respectively identified.

References
Shimadzu, H. (2018). On species richness and rarefaction: size- and coverage-based techniques
quantify different characteristics of richness change in biodiversity, Journal of Mathematical
Biology 77: 1363–1381. http://dx.doi.org/10.1007/s00285-018-1255-5.
Shimadzu, H., Dornelas, M. and Magurran, A. E. (2015). Measuring temporal turnoverin ecological communities, Methods in Ecology and Evolution 6(12): 1384–1394.
http://dx.doi.org/10.1111/2041-210x.12438.
2定量的生態学と統計モデル：
食われる側からみた特定食う種の貢献度
島谷健一郎 (統計数理研究所）
特定の天敵による死亡は、食われる側の死亡の中でどのくらいの割合を占めるか。Hoso and Shimatani (2020) では、野外における標識・再捕獲データと実験室で得た捕食成功データを階層ベイズモデルの中で統合することで、あるカタツムリ種が、天敵の１つであるヘビに、成熟して繁殖に至るまでどのくらい捕食されるかを定量的に推定した。さらに、ヘビに対してはトカゲのようにシッポを自切することで捕食を回避する場合もあり、それが成熟数をどのくらい増やす
かにより、進化的適応度を推定した。
このような定量的評価は、科学における営みの基本のはずである。にもかかわらず、現実にはこうした成果はなかなか伝統ある生態学の学術誌には掲載されない。その原因として、以下のような理由が考えられる。
１． 生態学が自然科学として未成熟で、多くの仮説や法則性は定性的記述にとどまる。定量的生態学の必要性は認識されていても、どう実践する
か、科学方法論（科学哲学）が確立されていない。
２． 統計モデルを用いて未知数量を推定したときの科学方法論（科学哲学）が確立していない。そのため、同様な問題は、統計モデルを用いた推定を要する科学全般に見られる。
３． 両者を融合させた、生態学において統計モデルを使って定量的仮説を検証する科学方法論（科学哲学）が確立されていない以前に、その必要性が認識されていない。
国内ではあまり議論されていないが、生態学の科学哲学は古くから議論されている。ただ、その多くは数理生態学に軸足を置いており、生態学における数理モデルの役割や、物理学と同じような科学方法論が適用可能か、などを論じている。2000 年代以降、先端の統計モデルは生態学で普通に適用されているが、そんな統計学を含めて論じている生態学の科学哲学論文は見当たらない。
統計学の科学哲学も古くから研究されている。ただその大半は、今日に至るも
なお古典検定と古典ベイズ（事前共役分布を用いる場合など事後確率が計算可能な場合）に関するもので、今日の階層ベイズモデルや状態空間モデルに関する
科学哲学は見当たらない。それどころか、最近のマルコフ連鎖モンテカルロ法

（MCMC）を用いるベイズモデルと古典ベイズを混同したような論文が公表されている。
逆に、統計学の業界から、階層ベイズモデルを用いる科学哲学も提起されている。ただ、そこでベースとされるのは、ポパー（反証可能な仮説が科学的仮説、などを主張）の時代、すなわちフィッシャーらにより統計学が確立されていく時
代以前の科学哲学にとどまりがちで、その後 50 年以上に及ぶ科学哲学の進展を踏まえた内容とは言い難い。
このような状況では、先端統計モデルを用いる生態学の科学哲学に関する議論が喚起されなくて当然である。このような現状についてまとめたのが、島谷（2021）や Shimatani (2021) である。
本講演では、最初に階層ベイズモデルを用いて定量的評価を行った研究事例を紹介し、それから、統計モデルを用いる生態学の科学哲学の未成熟さについて検討し、現在進行中の試みと展望について議論した。
参考文献
Hoso Masaki, Shimatani K. Ichiro. (2020) Life-History Modeling Reveals the Ecologicaland Evolutionary Significance of Autotomy. The American Naturalist 196: 690-703.
島谷健一郎 (2021) 推定目的の統計モデルを用いる帰納推論における数学・演
繹の位置付けと役割. Linkage.
Shimatani, I.K. (2021) Philosophy of statistical sciences: The roles of mathematics andstatistical model in estimation and other inductive inferences. 科学基礎論学会誌

ニワトリ雛の採餌⾏動解析―理論最適解から乖離する理由を探す―
統計数理研究所 川森愛[背景]
⾔語を持たない動物が何を考えて⾏動を決定するのか，その⽣々しい思考過程を知ることは容易ではない．しかし，思考の産物である⾏動を説明することによって，ある程度窺い知ることはできるはずである． 本発表ではニワトリ雛の採餌⾏動データを対象とし，複数の認知モデルを⽐較検討することで認知メカニズムの探究を試みる．
採餌における基本の⾏動原理は効率の最⼤化である．特に，餌場がパッチ状に分布する時の餌場滞在時間を考えるモデルとして，最適餌パッチ利⽤モデル(Marginal Value Theorem,
以下 MVT と表記)が提唱されている[1]．餌場にある餌は有限であるため，探索とともに発⾒に時間がかかるようになる．この状態を収益逓減と呼ぶ．⼀⽅，餌場間の移動中は収益 0の⼤損状態なので，頻繁に餌場を変えるのも得策ではない．そのため，⼀つの餌場にいつまで滞在するかの最適点が問題となるのである．
今，収益逓減のルールとして餌と餌の時間間隔が指数関数的に増加すると考える．具体的に， n 個⽬と n-1 個⽬の餌の間隔が arn-2 の指数関数で増加するとする．この時，MVT から予想される最適滞在時間 S*を求める式は⼤変複雑であり，このような計算を実際の動物が実装できているかは疑問である．しかしながら，⾏動データを⾒る限り理論を全く無視しているわけでもない．何らかの簡易的プロセスを⽤いて近似計算を⾏い，だいたいの最適性を実現していると考える⽅が現実的である．
[データ]
⾏動実験では，両端に餌場がある I 字型装置内を⽤い，ニワトリ雛に⾃由探索させた．餌場に到達すると⾃動的に粟が⼀粒ずつ供給されるが，供給までの待ち時間は指数関数的に増加する．装置の⾜場はトレッドミルとなっており，移動距離を⼈為的に 3 種類の⻑さに操作した．このように MVT が想定する餌環境を再現した実験で，ニワトリ雛が餌場に⼊ってから離脱するまでの滞在時間とその餌場に⼊るまでの移動時間の間の関係を調べた．ニ
ワトリ雛の装置内座標（1 次元データ）は超⾳波センサーにより⾃動的に，ニワトリ雛の⾏動とは無関係におよそ⼀定間隔で記録された．従ってデータとしては，滞在開始時刻として餌供給開始時刻(f)，餌場内での最終観測時刻(l)および餌場外での最初期観測時刻(m)の 3 つの時刻がある．また，離脱から次の滞在開始までの時間を移動時間(t)として算出した．これ
らのデータを⽤いてヒヨコの離脱時刻を説明する確率的意思決定モデルを作成し，⾏動データとのフィッティングを調べた．
[⾏動データへのフィッティング⽐較]作成したモデルに対しパラメータを最尤推定した上で，移動時間 1 点につき 100 回のシ

ミュレーションを⾏うことにより，滞在時間の予測値を求めた．最も成功したモデルの予測結果を⼀例として図に⽰す．⾚線が予測中央値，上下のピンク線が予測 95％範囲，⻘線は最適滞在時間の近似値，⿊丸は実際の⾏動データを⽰している．⾏動データは移動時間とともに分散が増⼤していく特徴がみられるが，このモデルではその様⼦をうまく捉えられて

滞在時間

いる．

予測中央値
平均利益率最⼤値

移動時間（図）モデルに基づいた滞在時間予測
複数の認知アルゴリズムを仮定したモデルを作成し検証した結果，ニワトリ雛が理論通りの最適滞在時間そのものを計算できている可能性は低いと判断できた．また，採餌効率を随時モニタリングするアルゴリズムにより滞在時間を決めると想定したモデルは全般的に，⾏動の分散を説明できる点で良いモデルであった．さらに，そのようなモデルの中でもどの変数を使⽤するかが重要であり，瞬間利益率のように認知的負荷が⼩さい指標を使⽤するモデルがニワトリ雛の⾏動にはよく当てはまった．
参考⽂献[1] E. L. Charnov (1976) “Optimal foraging: the marginal value theorem”, Theoretical Population
Biology, 9:129‒136

A uniﬁed formulation of k-means, fuzzy c-means andGaussian mixture model by the Kolmogorov-Nagumo
averageOsamu Komori
Seikei University
Kichijoji Kitamachi, Musashino, Tokyo 180-8633, Japankomori@st.seikei.ac.jp
Shinto Eguchi
The Institute of Statistical Mathematics
Midori-cho, Tachikawa Tokyo 190-8562, Japaneguchi@ism.ac.jp
October 26, 2021

Abstract
Clustering is a major unsupervised learning algorithm and is widely applied in data
mining and statistical data analyses. Typical examples include k-means, fuzzy c-means
and Gaussian mixture models, which are categorized into hard, soft and model-based
clusterings, respectively. We propose a new clustering, called Pareto clustering, based on
the Kolmogorov-Nagumo average, which is deﬁned by a survival function of the Pareto
distribution. The proposed algorithm incorporates all the aforementioned clusterings plus
maximum-entropy clustering. We introduce a probabilistic framework for the proposed
method, in which the underlying distribution to give consistency is discussed. We build
the minorize-maximization algorithm to estimate the parameters in Pareto clustering. We
compare the performance with existing methods in simulation studies and in benchmark
dataset analyses to demonstrate its highly practical utilities.

keywords: k-means; fuzzy c-means; Gaussian mixture model; maximum-entropy clustering;
Kolmogorov-Nagumo average; generalized energy function; Pareto distribution

1

1

Introduction

In data analysis or data mining, there are two fundamental types of methodologies: clusteringand classiﬁcation (Rokach & Maimon, 2005). Clustering, which is categorized as an exploratory
paradigm, detects the underlying structure behind the data and grasps the rough image before
proceeding to more intensive and comprehensive data analysis (Tukey, 1980; Dubes & A.K.Jain,
1980). On the other hand, classiﬁcation predicts unknown class labels of test data based onmodels constructed by training data with known class labels. The former is called supervised
learning, while the latter is called unsupervised learning in pattern recognition (Jain, 2010).
Clustering algorithms fall roughly into three categories: hierarchical, partitioning, and mixture model-based algorithms (Ghosh & Dubey, 2013). In hierarchical clustering, each observation is considered as one cluster in the initial setting. Then clusters are merged recursively
based on a similarity matrix deﬁned beforehand. The resultant clusters are expressed as adendrogram. The partition algorithm starts with a ﬁxed number of clusters and searches forthe cluster centers to minimize an objective function such as the squared distances between
the centers and observations. It ﬁnds the centers simultaneously. The model-based algorithm
assumes a mixture of probability distributions, which generates the observations, and assigns
the distributions to one of the mixture components. A Gaussian-mixture distribution-based
approach is widely used in this context.
In this paper, we propose a new clustering, called Pareto clustering in the framework ofquasi-linear modeling (Komori et al., 2016; Omae et al., 2017; Komori et al., 2017). It combines the cluster components by the Kolmogorov-Nagumo average (Naudts, 2011) in a ﬂexible
way. We consider a generalized energy function as an objective function to estimate cluster
parameters, which is an extension of the energy function proposed by Rose et al. (1990). The
objective function consists of a survival function of the Pareto distribution, which is widely usedin extreme value theory (Beirlant et al., 2004). We investigate the consistency of the parameters, resulting in the underlying probability distribution of the generalized energy function. We
ﬁnd that k-means (Cox, 1957; MacQueen, 1967) and fuzzy c-means (Bezdek et al., 1984) havethe underlying probability distributions with singular points at the cluster centers. This factshows a clear diﬀerence from the model-based clustering such as a Gaussian-mixture modeling. Moreover, we show that the quasi-linear modeling based on the the Kolmogorov-Nagumo
average connects k-means, fuzzy c-means and a Gaussian-mixture modeling using the hyper
parameters of the generalized energy function. See (Hathaway & Bezdek, 1995; Yu, 2005) for
the discussion of the relation between k-means and fuzzy c-means.

2