実践データ駆動科学オンラインセミナー
8th seminar データ駆動科学の最前線
複雑ネットワーク解析 / ベイズ統計

ベイズ推論に基づく物理モデル構築過程の効率化
統計数理研究所 統計的機械学習研究センター
特任助教 本武 陽一2021年7月30日
16:00 ~ 16:30
オンライン（ZOOM）

1

自己紹介
福岡県福岡市出身
学歴：
2008/3 東北大学

理学部

2010/3 北海道大学物理学科卒業大学院理学院

宇宙理学専攻修士課程修了 指導教官：石川健三 教授（素粒子物理学、テーマ：波束のダイナミクス）
2013/3 東京大学

大学院総合文化研究科 広域科学専攻

修士課程修了 指導教官：植田一博 教授（認知科学、テーマ：脳活動を用いたHAI評価）
2016/3 東京大学

大学院総合文化研究科 広域科学専攻

博士課程修了 指導教官：池上高志 教授（複雑系、テーマ：深層学習、群れ）
職歴：
2016/5 〜 2019/3 東京大学
特任研究員
2019/4 〜
特任助教

大学院新領域創成科学研究科

複雑理工学専攻（岡田真人教授、テーマ：データ駆動科学）
統計数理研究所

統計的機械学習研究センター（福水健次教授、テーマ：統計的機械学習）

2自己紹介現象3
武谷三段階論 [武谷 1942, 1966]

機械学習が得意

人間が得意

実体論

本質論

現象の内挿的モデル化と予測抽象化による大胆な外挿現象のモデル化

気体の分子運動論

一般原理の抽出

エントロピー増大則

dS = dQ/T

蒸気機関

機械学習⇄科学者
現象の予測
化学反応速度論

化学産業

3原理の転移化学ポテンシャルμ = ∂F/∂N

Willard Gibbs

自己紹介現象
大同特殊鋼 HPより

4
武谷三段階論 [武谷 1942, 1966]

機械学習が得意

人間が得意

実体論

本質論

現象の内挿的モデル化と予測抽象化による大胆な外挿現象のモデル化

気体の分子運動論

一般原理の抽出

エントロピー増大則

dS = dQ/T

磁性材料

機械学習⇄科学者

NEDO Project Success Story

構造材料

現象の予測
化学反応速度論

4原理の転移化学ポテンシャルμ = ∂F/∂N

Willard Gibbs

5

本日の内容
1. 科学者の営みとベイズ推論
2. ベイズ的線形回帰
3. ベイズ推論の利点
4. ベイズ推論の応用事例紹介
5. まとめ6

1. 科学者の営みとベイズ推論

1. 科学者の営みとベイズ推論

7

ベイズ統計
A

B
A∩B

条件付き確率●P(A | B) ==Thomas Bayesn(A ∩ B)n(B)n(A ∩ B)/N P(A ∩ B)=n(B)/NP(B)n(B ∩ A)
●P(B | A) =n(A)n(B ∩ A)/N P(B ∩ A)=
=n(A)/NP(A)

ベイズの定理P(A | B)P(B)P(B | A) =P(A)1. 科学者の営みとベイズ推論

8

ベイズ統計
因果則
ベイズの定理 ：P(A | B)

➡逆転 ：A → B
A: 結果
P(B)：原因の事前知識

A

原因
物理モデル・変数

結果測定・観測データ

事前知識 P(B)

因果関係 ：B → A
B: 原因

B

P(B | A)

B

A

原因
物理モデル・変数

結果
測定・観測データ

真実がゆらぐ！

P(A|B)：因果モデル
P(B|A)：原因の推定結果
データから物理モデル・変数の尤もらしさを推定！
真実はひとつではない1. 科学者の営みとベイズ推論

9バネの長さ: x

ベイズ統計
因果則

P(A | B)B

A

θ1 = {k, x0}

D = {xi, yi}Ni=1

M1 : y = k(x − x0)
M2 : y = k(x − x0)2

θ2 = {k, x0}

事前知識 P(B)

バネの初期長さx0

P(B | A)

B

A

θ1 = {k, x0}

D = {xi, yi}Ni=1

M1 : y = k(x − x0)
M2 : y = k(x − x0)2

力:y

θ2 = {k, x0}

バネの長さ: x

物理モデルや物理変数は一意に定まらない。
これが自然の摂理と考えるのは気持ち悪い！

1. 科学者の営みとベイズ推論

10

科学者（人間）の営みとベイズ推論
バネの長さ: x

観測・測定

科学者

事前知識

モデル候補
y = k(x − x0)

測定装置はあまり信頼できない

y = k(x − x0)2

実験者が初心者
信頼できる先生の発言
シンプルなモデルが良いモデル
エントロピー増大則
太陽は東から登る

バネの初期長さx0

地震が起こった
湿度が高かった
暑かった・寒かった

力: y

宇宙線が飛んできた
モデルのパラメータはだいたいこれくらいバネの長さ: x

y = 1.2(x − 1.12)「真実（犯人）は（80%の確率で）これだ！➡y = 1.2(x − 1.12)」

科学者の営みをベイズ推論としてモデル化すると考える11

2. フックの法則のベイズ推論

2. フックの法則のベイズ推論

12バネモデルの構築【知っていること】
因果則P(A | B)

A

M1 : y = k(x − x0)θ1 = {k, x0}

D = {xi, yi}Ni=1

M2 : y = k(x − x0)2
θ2 = {k, x0}(yi − k(xi − x0))2exp −2σ 2()
2πσ 21P(yi | k, x0, M1) =力:y

バネの長さ: x

因果則 P(Y | k, x0, Mk)

B

回帰式を中心に測定誤差（ガウスノイズ）が乗る

バネの長さ: x
バネの初期長さx0

P(A | B) = P(y1, y2, …, yN | k, x0, M1)
=

N

∏
i=1

=P(yi | k, x0, M1) =N21exp −( 2πσ 2 )

N

∏
i=1

(yi − k(xi − x0))2exp −2σ 2()
2πσ 21N∑i=1 (yi − k(xi − x0)2)
2σ 2

➡

最小二乗法

2. フックの法則のベイズ推論

13

ベイズ推論によるバネ定数推定【知りたいこと】
事後分布 P(k, x0 | Y)

「知ってること」：因果則（順モデル）

P(y1, y2, …, yN | k, x0, M1)
X

ベイズの定理 ：

M

k, x0

Y

因果関係 ：B → A➡逆転 ：A → B「知りたいこと」：尤もらしいバネ定数（パラメータ）
P(y1, y2, …, yN | k, x0, Mk)P(k, x0, Mk)P(k, x0 | y1, y2, …, yN, Mk) =P(y1, y2, …, yN )x0k

2. フックの法則のベイズ推論

14

ベイズ推論によるバネモデル評価【知りたいこと】
事後分布 P(Mk | Y)

「知ってること」：因果則（順モデル）

P(y1, y2, …, yN | k, x0, M1)
X

ベイズの定理 ：Mk, x0Y

因果関係 ：B → A

P(Mk | y1, y2, …, yN )➡逆転 ：A → B「知りたいこと」：モデルの尤もらしさ

P(y1, y2, …, yN | Mk)P(Mk)
P(Mk | y1, y2, …, yN ) =P(y1, y2, …, yN )M1

M2

P(y1, y2, …, yN | k, x0, Mk)P(k, x0 | Mk)P(Mk)= dkdx0∫
P(y1, y2, …, yN )

15

3. ベイズ推論の利点

3. ベイズ推論の利点

16

ベイズ推論の利点1
バネの長さ: x

データにない物理学者の知識を事前分布として導入できるP(x0)y = k(x − x0)
力: yx0 < 0ox バネの長さ: x0x0

3. ベイズ推論の利点

17

ベイズ推論の利点1
バネの長さ: x

データにない物理学者の知識を事前分布として導入できるP(x0)y = k(x − x0)
力: yx0 < 0ox バネの長さ: x0x0
➡少ないサンプル数で高い推定精度のモデル構築が実現できる➡事前知識を定量的に議論する基盤になる

3. ベイズ推論の利点

18

ベイズ推論の利点2事後分布を通してパラメータ推定精度や物理モデルの性質がわかる●物理モデルの理解

パラメータ2●推定精度

パラメータ1➡パラメータ1とパラメータ2は相関関係

3. ベイズ推論の利点

ベイズ推論の利点3
このデータは何次関数？→3次関数

y = 10x 3 − 12x19

3. ベイズ推論の利点

ベイズ推論の利点3●線形回帰

20f(x) = w0 + w1x + w2 x 2 + w3x 3 + w4 x 4 + w5x 5 + w6 x 6 + w7x 7 + w8x 8

3. ベイズ推論の利点

ベイズ推論の利点3●線形回帰

21f(x) = w0 + w1x + w2 x 2 + w3x 3 + w4 x 4 + w5x 5 + w6 x 6 + w7x 7 + w8x 8

3. ベイズ推論の利点

22

ベイズ推論の利点3
過適合を考慮したモデル評価ができる
最小二乗法（最尤法）➡二乗誤差を最小化
1 N2
E=y−f(x;w)( i)i2N ∑i=1w: モデルパラメータ

過適合を防ぐ−log P(Mk | y1, y2, ⋯, yN ) ∼Nyi − f(xi; w))(∑2i=1NE （エネルギ）※鞍点近似が成り立つ場合

d+log N2モデルの複雑さに対する罰則項

d: パラメータ数

3. ベイズ推論の利点

23

ベイズ推論の利点のまとめ
利点1：データにない科学者の知識を事前分布として導入可能

利点2：事後分布から物理量の推定精度や物理モデルの性質がわかる

利点3：過適合を考慮したモデル評価が可能

24

4. ベイズ推論の応用事例紹介

4. ベイズ推論の応用事例紹介

ベイズ推論による物理モデル構築の支援
同一の現象に対して複数の縮約モデル候補⇨ベイズ推論を用いて各モデルを評価・解釈
金属のクリープ過程の縮約モデル

NEDO Project Success Story[Mototake et al., Scientiﬁc Reports, 2020]

3d-X線光電子分光過程の縮約モデル

SPring-8 [SPring-8 HPより]

日本分析機器工業会 HPより[Mototake et al., JPSJ, 88(3), 034004, 2019]25

4. ベイズ推論の応用事例紹介

26

A universal Bayesian inference framework
for complicated creep constitutive equations
Yoh-ichi Mototake1, Hitoshi Izuno2, Kenji Nagata2,Masahiko Demura2, and Masato Okada3
1The Institute of Statistical Mathematics
2Research and Services Division of Materials Data and Integrated System, National Institute for Materials Science
3Graduate School of Frontier Sciences, University of Tokyo

Scientific Reports 10, 1, (2020).

4. ベイズ推論の応用事例紹介

27

金属のクリープ破壊●金属のクリープ現象●クリープ試験高温・高圧下での構造材料の破壊現象

http://fykai.blog.fc2.com/blog-entry-140.html

NEDO Project Success Story

Acquisition, Technology & Logistics Agency

http://fykai.blog.fc2.com/blog-entry-140.html
NIMS HPより
https://www.nims.go.jp/news/archive/2020/02/202002280.html

4. ベイズ推論の応用事例紹介

Scientific Reports 10, 1, (2020).
28

金属のクリープ破壊●クリープ試験0年30年

4. ベイズ推論の応用事例紹介

Scientific Reports 10, 1, (2020).
29

クリープ構成式によるアプローチ
複数提案されているクリープモデル
𝜀 = 𝜀0 + 𝑤1(1 − exp(−𝛼𝑡)) + 𝑤2𝑡 + 𝑤3(exp(𝛽𝑡) − 1)𝜀 = 𝜀0 + 𝑤1(1 − exp(−𝛼𝑡)) + 𝑤3(exp(𝛽𝑡) − 1)・・・

どれが適切なモデルなのか？

http://fykai.blog.fc2.com/blog-entry-140.html

http://fykai.blog.fc2.com/blog-entry-140.html

4. ベイズ推論の応用事例紹介

Scientific Reports 10, 1, (2020).
30

クリープ構成式によるアプローチ
1.クリープ曲線を構成式で最小二乗回帰

2.回帰パラメータの温度・応力依存性を回帰ϵkim(t; θkim) = ϵ0 + a1t b1 + a2t b2 + c1 exp(d1t) + c2 exp(d2t)最小二乗法 測定ノイズは考慮していない（そもそも計測系の測定ノイズは不明）
http://www.expertsmind.com/topic/material-problems-at-elevated-temperatures/effect-of-temperature-and-stress-919147.aspx[Kimura et al. 2014, 2015]
https://asmedigitalcollection.asme.org/PVP/proceedings-abstract/PVP2018/51678/V06AT06A059/278068
https://asmedigitalcollection.asme.org/PVP/proceedings-abstract/PVP2015/57007/V06BT06A058/285105

4. ベイズ推論の応用事例紹介

Scientific Reports 10, 1, (2020).
31

木村モデル
木村モデル：state of the artsなモデル
b1b2ϵkim(t; θkim) = ϵ0 + a1t + a2t + c1 exp(d1t) + c2 exp(d2t)[Kimura et al. 2014, 2015]

b2はほぼ1?
C2[Kimura et al. 2014]

4. ベイズ推論の応用事例紹介

Scientific Reports 10, 1, (2020).
32

木村モデル
木村モデル：state of the artsなモデル
b1b2ϵkim(t; θkim) = ϵ0 + a1t + a2t + c1 exp(d1t) + c2 exp(d2t)θkim = {a1, a2, b1, b2, c1, c2, d1, d2}[Kimura et al. 2014, 2015]

モデル比較
b2はほぼ1?

修正木村モデルb1ϵmod(t; θmod) = ϵ0 + a1t + a2t + c1 exp(d1t) + c2 exp(d2t)θmod = {a1, a2, b1, c1, c2, d1, d2}

4. ベイズ推論の応用事例紹介

Scientific Reports 10, 1, (2020).
33

提案モデルによる回帰結果
木村モデル

修正木村モデル

どちらも高精度の回帰結果かつ、わずかな差しかない➡測定誤差の値が重要になる

しかし、測定誤差はレンジでしかわからない➡測定誤差を確率変数としたベイズ推論

4. ベイズ推論の応用事例紹介

Scientific Reports 10, 1, (2020).
34

ベイズ推論に期待すること
利点1：データにない科学者の知識を事前分布として導入可能➡レンジで与えられた測定誤差の情報をモデル評価に取り込む

利点2：事後分布から物理量の推定精度や物理モデルの性質がわかる➡木村モデルでなぜC2の温度・応力依存性が複雑化するのか解明する

利点3：過適合を考慮したモデル評価が可能➡修正木村モデルが木村モデルより良い物理モデルか比較する

4. ベイズ推論の応用事例紹介

Scientific Reports 10, 1, (2020).
35

ベイズ的クリープモデル選択
b1b2ϵkim(t; θkim) = ϵ0 + a1t + a2t + c1 exp(d1t) + c2 exp(d2t)σ【知っていること】
因果則（順モデル）【知りたいこと】

θkimuの事後分布

ϵθkim

Mkt
N
21P(ϵ1, ϵ2, …, ϵN | θkim, σ, Mkimu) =exp −( 2πσ 2 )N∑i=1 [ϵi − ϵkim(ti; θkim)]22σ 2P(ϵ1, ϵ2, …, ϵN | θkimu, σ, Mkimu)P(θkimu | Mkimu)P(σ | Mkimu)
P(θkimu, σ | ϵ1, ϵ2, …, ϵN , Mkimu) =P(ϵ1, ϵ2, …, ϵN )【知りたいこと】
P(ϵ1, ϵ2, …, ϵN | θkimu, σ, Mk)P(θkimu | Mk)P(σ | Mk)P(Mk)
P(Mkimu | ϵ1, ϵ2, …, ϵN ) = dθkimudσ∫P(y1, y2, …, yN )モデルの事後分布

4. ベイズ推論の応用事例紹介

Scientific Reports 10, 1, (2020).

σの事前分布とモデル選択結果
ゲージ目盛りから見積もられたおおよその測定誤差レンジ

36σの事前分布によらず誤差は木村モデルが常に小さい➡誤差の大小ではないモデル比較が実現➡測定誤差の事前知識に応じて選択されるモデルが切り替わる➡モデルの正確な評価には測定誤差の正確な見積もりが必要4. ベイズ推論の応用事例紹介

Scientific Reports 10, 1, (2020).
37

事後分布の分析b1b2ϵkim(t; θkim) = ϵ0 + a1t + a2t + c1 exp(d1t) + c2 exp(d2t)木村モデル

修正木村モデル木村モデル

修正木村モデル➡事後分布は、木村モデルが過学習を行っていることを示唆4. ベイズ推論の応用事例紹介

Scientific Reports 10, 1, (2020).
38

C2の温度応力依存性が複雑な要因
b1b2ϵkim(t; θkim) = ϵ0 + a1t + a2t + c1 exp(d1t) + c2 exp(d2t)➡修正木村モデルが温度・応力方向の推定において、木村モデルを改善できそう4. ベイズ推論の応用事例紹介

Scientific Reports 10, 1, (2020).
39

C2の温度応力依存性が複雑な要因
b1b2ϵkim(t; θkim) = ϵ0 + a1t + a2t + c1 exp(d1t) + c2 exp(d2t)第2項の冪が0より大きくても，2次クリープ領域でほぼ線形になる➡加速クリープ領域を3つの基底で回帰した結果過学習している可能性

40

5. まとめ5. 本日のまとめ

41

まとめ• ベイズ推論は、科学者の知的活動を自然にモデル化する。
• ベイズ推論には、科学者にとって有用な特性がある。
• ベイズ推論は我々が現象や物理モデルを理解することを助けてくれるツールとなる。
一方で、、、P(ϵ1, ϵ2, …, ϵN | θkimu, σ, Mkimu)P(θkimu | Mkimu)P(σ | Mkimu)
P(θkimu, σ | ϵ1, ϵ2, …, ϵN , Mkimu) =（１）P(ϵ1, ϵ2, …, ϵN )P(ϵ1, ϵ2, …, ϵN | θkimu, σ, Mk )P(θkimu | Mk )P(σ | Mk )P(Mk )P(Mkimu | ϵ1, ϵ2, …, ϵN ) = dθkimudσ（２）∫P(y1, y2, …, yN )複雑なモデルの場合、(1)事後分布の関数形を調べたり、(2)周辺化の積分を実行するのに数値計算（モンテカルロサンプリング）が必要になる。
➡職人技的な技術が必要になる & 真実はサンプリングだけからは確定できない➡お気軽にご相談ください（@営業）mototake@ism.ac.jp