《講義ノート》

物性研究・電子版 Vol. 8, No. 1, 081201（2020年2月号）

熊本大学・パルスパワー科学研究所 赤井一郎

測定回数を増やすと測定値のヒストグラムは図 2 の様になり、ある分布関数に漸近する。
N =100, Step: 0.005 W
N =10000, Step: 0.001 W

要旨
データ駆動科学は、ベイズ推定やスパースモデリング等の機械学習の手法を計測・解析に組み込み、従来の知見を超える情報を抽出する取り組みである。これまで科学計測は、現象の源に原因があり、因果律に従って結果が決まると考える一方、誤差論では、計測値は偶然誤差によって確率分布すると考えてきた。データ駆動科学で用いるベイズ推定ではその考え方を逆転させ、計測で得られたデータを「起点」として原因が確率的に分布すると考える。
一方スパースモデリングは適切な物理モデルを用いて、データに含まれる要素を顕著な主成分から順に抽出し、適切な規範に沿って現象を理解する正則化の方法である。
データ駆動科学の基本的思想は、計測データを最善の方法 (機械学習) で骨の髄まで解析して虚心坦懐にデータに対峙することである。その結果、一切の作為的な先入観を排除するとともに、例え院生であっても経験豊富な研究者と対等に戦うことが可能となる。特にベイズ推定では物性パラメータの統計的分布の評価が可能で、新しく発見した現象や物性の、誰もが納得できる統計的証拠を示すことができる。更に、一切の先入観を排除してデータだけから現象を説明する物理モデルの選択も可能である。講義では、従来の誤差論を復習し、ベイズ推定とスパースモデリングの方法論の紹介と、具体的に物性研究に適用した例を示し、若手研究者の参画を喚起したい。

1

復習:誤差論1.1

真値、最確値、誤差

最確値とその誤差
複数回の測定等を行なった後、後述する方法を用いて真値にもっとも近いであろう値を求めることが可能である。それを最確値 (the most probable value) という。しかし最確値も誤差が伴う。よって最確値は最確値の誤差を伴って次の様に表す。
最確値の表し方 : 最確値 ± 最確値の精度(2)
偶然誤差と系統誤差
誤差には偶然誤差 (random error) と系統誤差 (systematic
error) がある。
• 偶然誤差 ξr :実験において避けがたいもので、複数回の測定に対し真値のまわりに正負同等にばらついて値がふれる誤差 (図 1-(a))。
• 系統誤差 ξs :
測定において真値から系統的にずれた値を与える誤差 (図 1-(b))。
誤差 ξ はこれらの和になる。
ξ = ξr + ξs
(3)

(b)

6
4
2
0

160
120
80
40
0

4.56 4.60 4.64 4.68
R ( W)

4.56 4.60 4.64 4.68
R ( W)(b)(a)

図 2: ヒストグラム。 (a)100 回測定、(b)10000 回測定。

その分布関数がガウス分布 (正規分布) である。
ガウス分布
真値 X である物理量の測定値の分布は、x という測定値を得る確率が][1(x − X)2(4)
f (x) = √ exp −
2σ 2
σ 2π
となるガウス分布 (Gaussian distribution) に従う。
σ は分布の広がりを表わし標準偏差と呼ばれる。
s
この関数は次の特徴を持つ。
• 真値 X に対して対称。
• 真値 X で最大値をとる。
X• |x − X| ≫ σ で急激に 0 に近づく。
図 3: ガウス分布• この関数は分布確率なので、測定値 x の取りえる範囲
x = −∞ ∼ +∞ の全積分で1 となる。
∫∞

誤差
測定値の誤差 (error) ξ は、測定値 x の真値 X からのずれとして定義される。
誤差 = 測定値 − 真値, (ξ = x − X)
(1)
真値は分からないので、そのままでは誤差も分からない。

(a)

8

1=

真値
ある物理量を測定したとき、得られた値がどこまで正しいか、その評価が重要である。測定には必ず誤差 (不確かさ) が伴うため、測定値は「真値」に一致するわけではない。

真値

ガウス分布 (正規分布)

Distribution

物性科学研究のためのデータ駆動科学入門

1.2

Distribution

第 64 回物性若手夏の学校・講義 2019/8/6∼8/8

第 1 図: 測定値のばらつきの様子。(a) 偶然誤差のみの場合。(b) 偶然誤差と系
統誤差がある場合。

• 次の関係が成り立つ。
∫
σ2 =

∞

−∞

1.3

f (x) dx(5)(x − X)2 f (x) dx

(6)

−∞

間接測定と誤差の伝播

直接測定 (direct measurement) できる物理量 x, y から、ある物理量 z がある関数 f (x, y) で、z = f (x, y) と求められるとする。これを物理量 z の物理量 x, y による間接測定(indirect measurement) という。
平均値を用いた間接測定値の導出
直接測定可能な物理量 a, b, · · · と z に
z = f (a, b, · · · )の関係があるとする。
• 物理量 a: 平均値 a, 真値 A, 測定値 a1 , a2 , · · · an
• 物理量 b: 平均値 b, 真値 B, 測定値 b1 , b2 , · · · bn
物理量 z の最確値 z は、直接測定可能な物理量の平均値を用いて、
最確値 : z = f (a, b, · · · )(7)で得ることが出来る。
では、得られた間接測定値の誤差をどう評価すればよいであろうか?
測定可能な物理量 x, y から、z
関数 f (x, y) で物理量 z が次式Zで得られるとする。
Z+DZ
z = f (x, y)
x, y の真値を X, Y とし、そz=f (x,y)Y
れらの値から求まる z の真値 Z
yはX
Dx
Z = f (X, Y )
x
Dy
とする。
x, y が、それらの真値からそれぞれ ∆x, ∆y だけ微小変化
した時、物理量 z の微小変化 ∆Z は、偏微分を用いて、∂ f (x, y)∂ f (x, y)
∆Z =
(8)
∆x +
∆y
∂x
∂y
Z
Z

《講義ノート》

物性研究・電子版 Vol. 8, No. 1, 081201（2020年2月号）

εzm = 0.6745σzm
√(
)2
)2
(
∂f
∂f
ε2am +
ε2bm + · · ·
=
∂a Z
∂b Z

i=1

i=1

最小二乗法は直線回帰でよく使われる。n 個の測定値の対 (x1 , y1 ), (x2 , y2 ), · · · , (xn , yn ) が得られたとし、これらの値において以下が満たされるとする。
• yi の各値の確率誤差が等しく、xi には誤差が無い。
これらの測定値の対に次の直線関係が期待されるとする。
y = ax + b
直線回帰W

y:R( )

0.2

0.0
0.0

di

0.1
x : length (m)

0.2

a0
b0
b

σx = 0 なので、yi の測定
値と理論値 (y = axi + b) の
残差 δi を取り扱う。
δi = yi − (axi + b) (12)
この残差の 2 乗 δi2 の和 Sn∑S(a, b) =
δi2
(13)
i=1

S

a

GlobalMinimum

Local
Minimaq初期値

q初期値 qGoal
パラメータ q = {a, b, ...}

が最小値となる条件から得られる a, b が最確値である。最小となるには、(14) 式の連立条件が満たされればよい。
∂S
∂S
= 0,
= 0 (14)∂a∂b
この方程式を正規方程式 (normal equation) という。

任意の関数 f (x; a, b, · · · ) への回帰は同様に以下のように行うが、推定したいパラメータが関数に非線形で含まれる場合 (一般的にそうだが)、図 4 に示した様に (16) 式の解がユニークに決まらない。

(b)

ベイズ推定

2

が最小となる必要がある。この様に誤差の 2 乗の和を最小にして、最確値を推定する方法が最小二乗法である。

(16)

図 4: 最小二乗法の困難さ

(10)

最小二乗法
測定値の確率分布 g(X) は X が本当の真値であるときに最大となり、f (ξ) は正しいガウス分布となる。g(X) が最大になるには、誤差の 2 乗の和 (指数関数の中)nn∑∑(xi − X)2
(11)
ξi2 =
S=

0.1

i=1

に対して (16) 式の正規方程式を解いて行われる。
∂S∂S== ··· = 0
∂a
∂b

(a)

最小二乗法

1.4

関数回帰
σx = 0、y 値の誤差は同じ (σy = σyi ̸= 0) とした場合、任意の関数 f (x; a, b, · · · ) への回帰は次式の Sn∑2S=[yi − f (xi ; a, b, · · · )](15)

平均二乗誤差 S(a,b,... )∂ f (x, y)∂ f (x, y),
は、f (x, y) の
∂x
∂y
Z
Z
x = X, y = Y における偏微分の微係数である。
平均値を用いた間接測定値の誤差 (誤差の伝播)
直接測定可能な物理量 a, b, · · · と、得たい物理量 z に、
z = f (a, b, · · · ) の関係があり、物理量 a, b, · · · が互いに独立であるとする。
• 物理量 a :平均値 a, 真値 A, 平均値の分布の標準偏差 σam• 物理量 b :平均値 b, 真値 B, 平均値の分布の標準偏差 σbm
物理量 z の最確値 z [(7) 式] の分布の標準偏差 σzm は(9) 式で、確率誤差 εzm は (10) 式で求められる。これを誤差の伝播 (Propagation of error) という。
√(
)2
(
)2
∂f
∂f
2 + ···
2
σzm =
σam +
σbm
(9)
∂a Z
∂b Z
である。ここで、

データ駆動科学は、これまでのデータ解析の価値観を一変させるもので、以下のことが可能となる。
• データを骨の髄まで解析する◦ 間接測定される物理量の統計的精度を評価◦ 物理モデルの妥当性の評価◦ データ解析における後ろめたさを一掃する◦ 既存の解析法を凌駕する新解析法• 実験の最適化◦ 極限測定で、必要な精度を得るための実験条件の推定◦ 実験を行う前の事前シミュレーションによる実験の効率化• 理論研究へのフィードバック◦ 解析に用いる物理モデルの良し悪しを判定
本節ではデータ駆動科学で用いられるベイズ推定の概要をまとめる。

2.1

ベイズの定理まず、ベイズの定理 [1] とベイズ推定の表現に用いられる、同時確率 (joint probability) と条件付き確率 (conditional
probability) について述べる。
同時確率とは、A という事象と B という事象が同時に起こる確率 [図 5-(a)] で、P (A, B) ≡ P (A ∩ B) と表す。
一方、A という事象と B という事象が同時に起こる確率は、まず A という事象が起こる確率 P (A) と、A という事
象が起こった条件下で B という事象が起こる条件付き確率の積 P (A, B) = P (B|A)P (A) で書くことができる。ここで P (B|A) が、A という事象が起こった条件下で B という事象が起こる条件付き確率で、次の様に書ける。
U

A

U

A

U

A

P(A,B) B

P(A,B) B

P(A,B) B

(a)

(b)

(c)

図 5: 同時確率 (a) と条件付き確率 (b), (c)•A 条件下の B の確率 [図 5-(b)],•B 条件下の A の確率 [図 5-(c)],P (A, B)
P (A)
P (A, B)
P (A|B) =P (B)P (B|A) =《講義ノート》

物性研究・電子版 Vol. 8, No. 1, 081201（2020年2月号）

ベイズの定理布 (ガウス分布) に従う場合について考える。
この場合、i 番目の測定点で正規分布ノイズ ni が発生すU
200 年以上前に T. Bayes は、[1] で、 る確率 P (n ) は、正規分布により、
Ai()
同時確率について、以下のベイズのP(A,B) Bn2i1定理を提唱した。
√ exp − 2
(22)
P (ni ) =2σdataσdata 2π
P (A, B) = P (B|A)P (A)である。ベイズ推定ではこの仮定に基づき因果律を遡る。
図 6: ベイズの定理と= P (A|B)P (B) (17)
データセットの横軸 xi は確定していて与えられていると同時確率しよう。そうすると、(21) 式より、ノイズ ni は、
P (B|A)
P (A|B)
P (A|B) =P (A), P (B|A) =P (B)ni = yi − f (xi ; θ)
P (B)
P (A)
であるから、(22) 式のノイズ ni を与える確率 P (ni ) は、
(18)(20) 式の因果律に従って、パラメータセット θ が決まった(18) 式で表したベイズの定理について再認識をしておく。
条件下でデータ Di が与えられる条件付き確率
ベイズの定理の言うことは、(19) 式に示すように、事前条
P (ni ) = P (Di |θ)
件がない状況では A 事象が P (A) の確率で起こるのに対し、と等価であることが分かる。よって(21) 式を用いて ni を残その確率 P (A) が、B が起こった条件下では A の起こる確
差
y
−
f
(x
;
θ)
で書き換えると、(23)
式と書ける。
ii
率が、尤度 P (B|A)/P (B) によって、P (A) から P (A|B) へ
P
(n
)
=
P
(D
|θ)
i
i
と更新されることを意味する。
)(
2
事前確率、事後確率と尤度
1[yi − f (xi ; θ)]√ exp −=
(23)
2
P (B|A)
2σdataσdata 2π
P (A|B) =
P (A)
(19)
| {z }
P (B) | {z }データに重畳するノイズ ni は i に関して独立であるか| {z } 事前確率
事後確率
ら、パラメータセット θ が与えられた条件下でデータセッ尤度
ト D が得られる条件付き確率 P (D|θ) は、(24) 式の平均二• 事前確率 (prior probability)
乗誤差
E(θ) を用いて、(25) 式となる。
P (A) は、何の条件もないときに A の事象が起こる確
N率で、事象 A の事前確率 (prior probability) と言う。
1 ∑
2
E(θ) =[yi − f (xi ; θ)](24)• 事後確率 (posterior probability)
2N i=1
一方 P (A|B) は、事象 B が起こった後の確率を意味し、
)N
(
)
(
N
1
事後確率 (posterior probability) と呼ばれる。
√
exp − 2 E(θ)(25)
P (D|θ) =• 尤度 (ゆうど: likelihood)σdata
σdata 2π
この事前確率と事後確率を結ぶ係数にあたる
P( θ , D )= P( θ ∣D ) P ( D)P (B|A)/P (B) を尤度 (ゆうど: likelihood) と呼ぶ。
測定行為は因果律に因果律を遡る従い、(20) 式が示した様
2.2 因果律と物理モデル
実験データセットDに、パラメータセット θ
基本方程式
因果律 (causality) は、結果には何らかの原因があり、そyが確率分布 P (θ) で決ま
因(物理モデル)の間に一定の法則性があることを言う。実験の測定行為に果った上でP (D|θ) を得る
パラメータ
律
おける因果律についてまとめておく。
セット qことである。その結果、
x
因果律と P (θ, D)P( θ , D )= P( D∣θ ) P ( θ )P (θ, D) の同時確率が得ある物理現象をられることに対応する。
実験データセットD
図 8: ベイズの定理の適用
測定することを考
基本方程式
y
一方ベイズの定理によれば、P (θ, D) は (26) 式と書ける。
えよう。
(物理モデル) 因
因果律を遡る :P (θ, D) = P (θ|D)P (D)(26)
物理現象にはそ果
パラメータ
律
れが従う基本方程
セット q
ベイズの定理の適用
x式があり、基本方
我々人類が得たいものはデータセット D が得られた条程式に従って測定
件下でのパラメータセット θ の確率分布であって、その量を与える物理モ
図 7: 測定行為の因果律
P (θ|D) はベイズの定理を用いて、(27) 式と書ける。
デルが存在する。
P (D|θ)物理モデルには、現象や測定量を特徴づける複数のパラ
P (θ|D) =
P (θ)
(27)
| {z }
P (D) | {z }メータ (パラメータセット θ) が内在する。
|
{z
}
事後確率
事前確率
因果律では、基本方程式、物理モデルとパラメータセッ
尤度ト θ が現象の原因であり、測定することによって得られ
尤度の P (D) はデータが取得された時点で値が確定する
るデータセット D が結果となる。
ため、P (θ|D) は (25) 式より、(28) 式となる。
()
物理モデルのパラメータセット θ に基づき、あるデー
N
P
(θ|D)
∝
exp
−
E(θ)
P (θ)
(28)
タセット D が得られる同時確率を考えよう。因果律に従2| {z }| {z }
σdata
えば、原因であるパラメータセット θ が確率分布 P (θ) で
事後確率
事前確率
決まっていて、測定することにより、P (D|θ) が得られ、2.4 データ駆動科学の最初の一歩次式のように同時確率が得られることになる。
表 1 にまとめた様に、誤差論とベイズ推定は根本的なデー
因果律 :P (θ, D) = P (D|θ) P (θ)(20)| {z }
タに対する認識が異なる。
←−−−−−
測定行為

2.3

データ駆動科学におけるベイズ推定[2] を参考文献として、話をベイズ推定に戻そう。データ
セット D は (21) 式の様に、
yi = f (xi ; θ) + ni
(21)
| {z } |{z}物理モデル

表 1: これまでの誤差論とベイズ推定の考え方の違い

f (xi ; θ)

θ

yi

誤差論

確定的

確率的

ベイズ推定

確率的

確定的

ノイズ

物理モデル f (xi ; θ) に偶然誤差 ni が重畳して与えられる。
以降このノイズが標準偏差 σdata で特徴づけられる正規分 • これまでの誤差論: 「パラメータ θ は真値で確定されて
いて、データがノイズで確率的に与えられる」と考える。

《講義ノート》

物性研究・電子版 Vol. 8, No. 1, 081201（2020年2月号）

マルコフ連鎖モンテカルロ法(28) 式の P (θ|D) をどの様にして得るか? ベイズの定理は200 年以上前に提唱されていたにもかかわらず、その実質的な計算には近年のコンピュータの進化を必要とした。
一番単純に P (θ|D) を評価するアルゴリズムは、パラメー
タ空間を必要な精度でメッシュに切って総当たりするグリッドサーチである。まずその例として、表 2 に示した様にパラメータは 5 つで、指定した範囲内を必要な有効桁の精度
で、P (θ|D) を評価することを考えよう。この場合 P (θ|D)
物理量

真値

範囲

有効桁

単位

f

3.0

2.5 ∼ 3.5

3桁

THz

ϕ

0.0

−π/2∼+π/2

3桁

rad.

τr

0.4

0.1 ∼ 1.0

2桁

ps

τd

6.0

1.0 ∼ 10.0

2桁ps

メトロポリス法マルコフ連鎖モンテカルロ法の代表的なアルゴリズムに
メトロポリス法 (Metropolis algorithm)[4] がある。実質はコンピュータの高速化により近年になって具現化 [5] された。
メトロポリス法は、以下の過程で P (θ|D) の確率分布を得るアルゴリズムである。
1. 乱数 ε を用いて、確率的に遷移候補を生成する。
θ′ = θt + ε,(31)
2. (28) 式から求められる採択確率 α
P (θ′ |D)
α = min{1, r}, r ≡(32)
P (θt |D)
で、θt+1 ← θ′ の状態更新を行う。採択されなかった場
合は、θt+1 ← θt 更新を行う。
特に候補値生成用の乱数 ε に、左右対称な一様乱数、正規乱数なども用いるものをメトロポリス法と呼ぶ。
その制限を拡張したものは、メトロポリス・ヘイスティン
グス法 (M-H 法: Metropolis-Hastings algorithm) と呼ばれる。
メトロポリス法の概要を以下に示す。
E(θ)

P(θ'|D)a
1.0
0.1 ∼ 2.0
2桁
表 2: 解析により得たい物理量とその有効桁
0

θ
θ

θ0

1

2

モンテカルロ法は、カジノで有名な国家モナコ公国にあるモンテカ
ルロの地名に由来する。カジノでサイコロを振るように、乱数を用いる方法であるためである。

0

θ
θ

θ0
θt

θt + εIterations

Iterations

θ

θ

θ'

t+1

θt+1

θ

θ': be accepted

+1

θ

(a)

(b)

図 9: メトロポリス法の状態更新の例-1

E(θ)

E(θ)

?
P(θt+1|D)?
P(θt+1|D)

1r

P(θ'|D)

P(θ'|D)P(θ|D)
0

θ
θ

θ0

t

θt

t+1

θt+1 + εθt+1P(θ|D)
0

θ

t

θ

t+1

θ
θ

θ0
θt
θt+1 + εθ'(c)
マルコフ連鎖は、A. A. Markov の「マルコフ過程」名に由来したものである。マルコフ過程は、物理現象を理解する確率過程で、未来の挙動が現在の挙動のみによって決定され、それ以前には依存しない確率過程のことである。

P(θ|D)t

Histogram

マルコフ連鎖モンテカルロ法 1 2 (MCMC: Markov chain
Monte Carlo methods)[3] が開発されている。ここではその基本的な考え方をまとめておく。
1. まず、各パラメータに「適当」な初期値 θinit. を与える。
真値 (不明) らしき値でも構わないし、乱数で適当にパラケさせても良い。
2. 次に、現在のパラメータの値 θt を元にして、何らかの乱数 ϵ によって、パラメータの次の候補 θ′ を生成する。
θ ′ = θt + ϵ
3. 現在の値 θt と、候補の値 θ′ の事後確率の比 r を評価する。
P (θ′ |D)r≡
(29)
P (θt |D)
ここで比を評価する意味は、(27) 式に現れる P (D) を評価する必要がなくなり、事前確率 P (θ) が一様で、 重畳するノイズが正規分布ノイズならば[]N′
r = exp − 2 {E(θ ) − E(θt )}
(30)
σdata
となることである。
4. 比 r の大小で、候補の値 θ′ の採択・非採択を決定し、次のパラメータ値 (2. から繰り返す) とする。
マルコフ連鎖モンテカルロ法では、2.∼4. の繰り返し過程において、マルコフ連鎖的にθinit → · · · → θt−1 → θt → θt+1 → · · ·P (θ|D) を与える度数分布を、繰り返しながら蓄積する。
よってマルコフ連鎖モンテカルロ法では、P (θ|D) を与える度数分布が得られた時点で繰り返しを終了し、得られた度数分布の度数 N (θ) から、P (θ|D) を得る。

P(θt|D)P(θ|D)

Iterationsの計算は 4.835 × 10 回必要で、例えば P (θ|D) の 1 回の計算に 16.4 ms(実測例) かかったとすると、全ての計算をし終えるのに、25,000 年の時間を費やす。
この計算時間爆発を解消するアルゴリズムとして様々な
13

E(θ)

Histogram

2.5

2.6

Iterations• ベイズ推定 : 「データが確定されていて、θ が P (θ|D)
によって確率的に決まる」と考える。
この様な「データを出発点と考える」考え方は、データ駆動科学の最初の一歩である。

θt+1

θ
θ

θ'

θ': be acceptedstochastically with r(d-0)

図 10: メトロポリス法の状態更新の例-2• 図 9-(a) : θ の初期値 θ0 を適当に設定する。
• 図 9-(b) :1. 現在の状態 θt から、乱数 ε で候補 θ′ (= θt + ε) を生成する。ε は、平均値 = 0 の一様乱数や正規乱数等。

《講義ノート》

物性研究・電子版 Vol. 8, No. 1, 081201（2020年2月号）
E(θ)

E(θ)
P(θt+1|D)1r

P(θt+1|D)

1Ur

P(θ'|D)

0

θ0

t

θt

t+1

θt+1 + εt+2Iterations

P(θ'|D)P(θ|D)

θt+1

P(θ|D)

θ
θ

0

θ

t

θ

t+1

θt+2

θt
θt+1 + εθ': be accepted

θ

θ
θ

θ0

θ'
t+2

Iterations

U

θt+1

θ

記号

θ̄ ± ϵθ振動数

f2930.04 ± 0.93 GHz

位相ϕ5.219 ± 0.012 rad.

立上り時定数τr

0.0809 ± 0.0048 ps

減衰時定数

τd

2.811 ± 0.051 psa(8.63 ± 0.10) × 10−16
表 3: パラメータセットとその平均値/Table I in [6].
積分強度

θ

θ'
θt+2

物理量

θ

θ': be rejected

る正規分布ノイズの標準偏差 σdata は、信号強度の次元で、
σdata = 3.1 × 10−5 と設定した。
+1
メトロポリス法の繰り返し回数は 100 万回である。その+1
内、初期の 10 万回は、適当に設定した初期値から適切なパラメータ値付近に収斂する焼きなまし区間 (Burn-In phase)θθとして破棄し、それ以降の度数分布を事後確率分布 P (θ|D)
(d-1)
(d-2)
として取得した。表 3 の ±ϵθ は P (θ|D) の分布幅を公差誤
差(確率誤差) で示したものである。
図 11: メトロポリス法の状態更新の例-3
このベイズ推定で注目すべき点は、周波数 f の推定精度である。従来、減衰振動波形の振動モードの抽出にはフーリエ変換等が用いられるが、振動振幅が減衰特性を持つた
2. 事後確率を比較する。
め、フーリエ変換スペクトルで不確定性幅は免れない。そ′P (θ |D)r=>1
れに対しこのベイズ推定の結果では、振動振幅の変化を適
P (θt |D)
切にモデル化することによって、フーリエ幅と比較して 2
図 9-(b) では、候補 θ′ の方が確率が高い (r > 1) ので、
桁程度高精度な周波数の推定に成功 [6] している。
無条件で候補を採用 (θt+1 ← θ′ ) して記録。
• 図 10-(c): 現在の状態 θt+1 から、乱数 ε で候補 θ′ (θ′ = 2.8 候補値を生成する乱数の大きさ σθθt+1 + ε) を生成する。
MCMC 法では、パラメータ候補値の生成で用いる乱数の• 図 10-(d-0): 事後確率を比較する。
大きさ σθ を適切に選択する必要がある。
P (θ′ |D)
原理的には得られるパラメータの P (θ|D) の広がり (幅)r=<1
は、実験データに重畳するノイズ ni の標準偏差 σdata によっ
P (θt+1 |D)
図 10-(d-0) では、r < 1 で、候補 θ′ の方が確率が低い。 て決定づけられるが、不適切な σθ 、例えば小さすぎる σθ• 図 11-(d-1): 0 ∼ 1 の一様乱数 U (0, 1) を発生させ、r と を用いると値の遷移幅が小さく、広いパラメータ空間内で十分な履歴を得る効率が低下する。一方大きすぎる σθ を用比較する。
図 11-(d-1) では、U < r であるので、候補 θ′ を採用 いると、パラメータ空間の全域探索には有効であるが、確率分布比 r による採択効率が落ち、やはり、十分な履歴を(θt+2 ← θ′ )。
• 図 11-(d-2): 0 ∼ 1 で発生させた一様乱数 U (0, 1) が、U > 得る効率が低下する。よって、メトロポリス法では、確率r であるので、候補 θ′ を棄却し、θt+1 を採用 ( θt+2 ← 分布比 r をモニターして、適切な σθ を選択するチューニン
グが必要である。
θt+1 ) する。
• 以上を繰り返し、θt の履歴のヒストグラムとして P (θ|D)
2.9 最小 2 乗法とメトロポリス法の比較
を得る。

-4

D

CP signal ( R/R) ×10

Steepest descent method

Metropolis methodE(θ)

メトロポリス法を用いたベイズ推定の例[6] に示したベイズ推定結果の概要をまとめる。
図 12 は、ビスマスのコ
ヒーレントフォノン信号のBi thinfilm
2例である。
こ れ は 、フェム ト 秒 の0
レーザーパルスを試料に照射させると、そのレーザー光によって、試料内の格子-2振動や分子振動が一斉に励振され、その振動に応じて-4
物質の反射率や透過率が周02468 期的に変化する現象を時間Time Delay (ps)
差のつけた弱いプロープ光図 12: ビスマス薄膜のコヒーレン
パルスで計測して得られる
トフォノン信号信号である。
このコヒーレントフォノン計測は、試料内の格子振動や

E(θ)

2.7

Uniform
randomnumberHistogram

Histogram

Uniform
randomnumberLocal
Min.

Global
Min.

θ

Local
Min.

GlobalMin.

θ

図 13: 最急降下法等による最小 2 乗法とメトロポリス法の比較

通常の最小二乗法は、(33) 式に示した、測定データ yi とフィッティング関数 g(xi ; θ) の平均二乗誤差 E(θ)
n1∑[yi − g(xi ; θ)]2
(33)
E(θ) =ni=1を最小化して、パラメータセット θ を決定する。また、最分子振動のダイナミクスを時系列変化として得られる方法 小化させる方法として最急降下法が用いられる場合が多い。
しかし、E(θ) がパラメータセットの θ 空間で図 13 に示したである。このデータに
[
(
)
(
)]
a
t
t
ように、多谷構造を取る場合、パラメータセットの初期値
Y (t) =exp −− exp −sin (2πf t + ϕ)τd − τrτdτr
によっては、多谷構造の内、大域的極小値 (global minimum)の物理モデルを適用し、メトロポリス法でベイズ推定した。 に到達できず、局所的極小値 (local minimum) にトラップさパラメータは表 3 に示した 5 個である。データに重畳す れることがある。

《講義ノート》

物性研究・電子版 Vol. 8, No. 1, 081201（2020年2月号）

それに対しメトロポリス法では、以下のメリットがある。
• 局所的極小値から解放される。
パラメターの候補を乱数生成するため、局所的極小値に局在せず、原理的には大域的極小値の探索が可能である。
• パラメータセット全空間の事後確率が得られる。
詳細つりあいが保証されるため、無限に時間をかければ、パラメータセットの全空間にわたってパラメータセットの事後確率 P (θ|D) を得ることができる。

• パラメータセットの確率分布を得ることができる
確率的にパラメータ候補を生成し、その度数分布から事後確率 P (θ|D) を得るため、大域的極小値近傍で、図 13に示したように、パラメータセットの分布を生成することが可能で、事前確率 P (θ) が一様ならば、P (θ|D) は)(N
P (θ|D) ∝ exp − 2 E(θ)σdataとなる。ここで σdata は、データに重畳するノイズ (正規分布ノイズ) の標準偏差で、
2
σdata∼ kB Tnoise(34)と考えれば、ボルツマン分布になっていることが理解で
きる。

2.10

レプリカ交換モンテカルロ法概論
平均二乗誤差 E(q )メトロポリス法では原理
Local
Global的に一意の確率分布を得る
Minimum
Minimum
ことが可能であるが、E(θ)の局所的極小値を脱して大sdata
域的極小値に達する (図 14)には時間を必要とする。
度数
その問題点を回避するア
分布
パラメータ qルゴリズムとして、レプリカ
交換モンテカルロ法 [7, 8, 2]
図 14: メトロポリス法の問題点(RXMC: Replica eXchange
Monte Carlo method) が開発された。
この拡張は、(34) 式に示した様に、E(θ) をエネルギーと2
捉え、σdataと Tnoise が対応しボルツマン分布と捉えること、イジングスピンにおける熱アニール効果をうまく組み合わせた着想 [8] を起点とする。
レプリカ交換モンテカルロ法
この考え方を発展させ、人為的に逆温度()β を導入する。
N
P (θ|D, β) ∝ exp − 2 βE(θ) P (θ)
(35)
σdata
• β = 1 : (28) 式と等価。
• 0 < β < 1 : 有効温度を上げたことに対応する。
であるが、レプリカ交換モン
テカルロ法では、0 < β ≤ 1で複数の β のレプリカ (図 15)を用意して大域的極小値を高速に探索する。
図 15: 異なる β のレプリカ
その探索は、• 各レプリカで、メトロポリス法で P (θ|D, β) を取得。
• 隣り合うレプリカ間で状態交換を確率的に行う。
を繰り返すことで実行する。

2.11

ベイズ自由エネルギーに基づくモデル選択

モデル選択とは、実験データを説明する物理モデルが複数存在す
る際に、どの物理モデルが統計的に妥当であるかを判断することである。以下の例がそれに当たる。

yx
図 16: 直線回帰？

• 直線回帰単純な例として、図 16 のようなデータに対し、
y = b, y = ax, y = ax + b, · · ·どの物理モデルを選択すべきかの判断は難しい。
• スペクトル分解 [2]計測された様々なスペクトルを成分分解する場合を考えよう。その計測データにピーク形状のスペクトル要素がいくつ含まれているか? 多分このような問いかけをしたら、要素数の数は既知であると答えるかもしれない。
しかし既知と思った要素数は正しいであろうか? 新しい
ピーク構造の発見を見落としていないか?
• ハミルトニアン選択 [9]計測されたデータの理解は様々な論争を引き起こす。根底にあるものはハミルトニアンで、データを説明するためには、ハミルトニアンのどの項が重要なのか? 複数のモデルハミルトニアンでどれが、データをよく説明するか?
ベイズ自由エネルギー最小化を情報量規準としたモデル
選択 [2, 9] は、これらの問題や問いかけに対し、データを起点としたデータ駆動科学的な解決策を提供する。
このモデル選択を念頭に、因果律を再考しよう。これまでは物理モデルが決まっているとして、それに内包するパラメータセット θ が因果律の原因と考えてきた。
因果律
しかし実際は、図 17 に示した様に、その前に物理モデルの選択が、
K
Dq
大元の原因であり、モデル選択さベイズ推定れてから、その物理モデルのパラ
メータセット θ によってデータセッ
図 17: モデル選択
ト D が与えられるのである。
ベイズ推定のモデル選択は、データセットが与えられた条件下で物理モデルである K の確率分布 P (K|D) を評価することである。
図 17 の因果律に沿って K, θ, D の同時確率 P (K, θ, D) を表すと、(36) 式となる。
P (K, θ, D) = P (D|θ, K)P (K)(36)
ここで、P (D|θ, K) は K と θ が決まった条件下でデータ
セット D が与えられる条件付き確率であるから (28) 式と等価で、(37) 式である。
(
)
N
∏
N
P (D|θ, K) =P (Di |θ, K) ∝ exp − 2 E(θ, K)
σdata
i=1
(37)
モデ ル選択 では P (K|D) を評 価し たいわ けだ か ら 、P (K, D) の同時確率について考えよう。この P (K, D) は、
P (K, θ, D) で θ を消去すれば良いので、P (K, θ, D) を全パ
ラメータ空間 θ に渡って積分すれば得られる。この様に特定の変数を消去することを周辺化(marginalization) という。
∫

P (K, D) = P (K, θ, D) dθ(38)
ここで β < 1 では有効温度を上げたことに意味する。そのような高温のレプリカ内でメトロポリス法を用いて探索す その結果ベイズの定理を適用し、P (K|D) は (39) 式となる。
P (D|K)る際、適切な採択率でサンプリングの効率化を行うため、P (K|D) =
P (K)
P (D)
候補値生成で用いられる乱数幅 (ステップ幅) を広げて探索()∫
する。よって高温では、図 15 に示した様に、E(θ) の局所構
N
∝ P (K) exp − 2 E(θ, K) P (θ) dθ
造を無視して、パラメータ空間を大域的に探索が行われる。
σdata
この様に β < 1 の高温で大域的に探索がなされるため、(39)
隣り合うレプリカ間で状態交換を確率的に行うことで、例
熱統計力学において分配関数
Z
は、逆温度
β
と状態を規
え低温 β = 1 で局所的極小値にトラップされても、一旦高∑温に状態交換がなされ、局所的極小値からの脱出が効率化 定するパラメータ θ を用いて、Z = θ exp (−βE(θ)) と定義される。これと (39) 式との対応を見ると、平均二乗誤差される。

《講義ノート》

物性研究・電子版 Vol. 8, No. 1, 081201（2020年2月号）

2と捉えると、
E(θ, K) をエネルギーと捉え、β ← N/σdata
情報科学の分配関数
Z()∫
N
Z ≡ exp − 2 E(θ, K) P (θ) dθ
(40)
σdata
を定義することは極めて有用である。
またこの分配関数 Z を用いてヘルムホルツの自由エネル
ギー F に対応するベイズ自由エネルギー [2] が、次式で定義される。 [∫()]
NF ∝ − lnexp − 2 E(θ, K) P (θ) dθ
(41)
σdata
モデル選択は熱統計力学の考え方に則り、ベイズ自由エ
ネルギー最小化の情報量規準で行うことができる。
更にベイズ推定では、モデルが選択される事後確率
P (K|D) が評価 [(39) 式] できることにも意味がある。単純にモデルの妥当性を 0/1 のデジタルで判別するのではなく、そのモデル選択の妥当性が定量的に評価できるのである。

3

スパースモデリング

この節では、ベイズ推定と異なる発想に基づくスパース
モデリング (SpM)[10] を紹介する。
何かのグループで取った記念写真を見た時を想像してほしい。その時、背景に何が写っているかを気にせずに誰が写っているか? 何人写っているか? 誰かの影に隠れているのは誰だ? という様なデータ処理が頭の中でなされている。
SpM は、人間が物体や現象を認知・識別する際に頭の中で行われる処理行為を、実験データの解析に適用するものである。つまり、実験データを与える物理現象は、より簡潔な物理モデルで説明可能であるという前提に立ち、簡潔に説明できる基底関数を選択して、その基底関数で表させる
少数の主要成分 (スパース解) の抽出を可能とする。

3.1

数学的準備

3.1.1 p 次ノルム
ノルム (norm) は、平面や空間における幾何学的ベクト
ルの「長さ」の概念の一般化で、⃗
x ベクトルの p 次ノルム
(1 ≤ p < ∞) は以下で定義される。
( n
)1/p
∑
√
p
∥⃗x∥p ≡|xi |= p |x1 |p + · · · + |xn |p

誤差を最小化した θ⃗ を求める行為である。よって、次式と
書ける。
(
)
2
ˆ
⃗
⃗
最小二乗法 θ = arg min⃗y − f (⃗x; θ)
2

⃗
θ

3.2

連立 1 次方程式を解く

植村誠先生の「天文屋のための HOW TO スパースモデリング」の第 1 章 [11] を参考に進める。
計測データから主成分を抽出する行為は、観測値から推定パラメータを求める逆問題である。SpM では計測データ
を M 個の要素を持つ ⃗y とし、抽出したいパラメータ ω⃗ がN 個の要素を「持ちうる」とし、その間を−→⃗y = X ω
⃗
(42)

 


y1
X11 X12 · · · X1N
ω1

 


 y2   X21 X22 · · · X2N   ω2 

 


 . = .
 .  (43).. . .
.. 
 ..   ..
.
. .  . 
.

 


yM

XM 1 XM 2 · · · XM N

ωN−→の線形写像で書き表す必要がある。ここで X は物理モデル⃗ を行列で表したものであり、既知であれば、この線
f (⃗x; θ)形写像は ω
⃗ に関する連立 1 次方程式であることが分かる。
• 最小二乗法:
観測点数がパラメータ数よりも多い場合 (M > N ) なら
ば、次式のように最小二乗法を使えば良い。
(
)
−
→ 2
ˆ
最小二乗法 ω⃗ = arg min ⃗y − X ω
⃗
ω
⃗

2

• 不良設定問題:
観測点数がパラメータ数よりも少ない場合 (M < N ) は
不良設定問題 (ill-posed problem) と呼ばれる。この場合
最小二乗法だけでは、解が多数存在し、その中から 1 つの解を選ぶためには他にも条件が必要になる。SpM の 1種である「L1 正則化」はその様な手法の 1 つである。

3.3

L1 正則化

前述の様に不良設定問題を考える。

3.3.1 0 次ノルム最小化
0 次ノルム最小化
i=1
ω
⃗ がスパースであると仮定できるとき、単純に考えるとSpM では、次の 1 次, 2 次ノルムと 0 次ノルムが問題となる。
−→⃗y = X ω⃗• 1 次ノルムの条件を満たす複数の ω⃗ の解の中から、非ゼロ要素の「個∥⃗x∥1 ≡ |x1 | + · · · + |xn |数」(0 次ノルム: ∥⃗ω ∥0 ) が最小になる ω⃗ を探すのが良さそ• 2 次ノルム√
うに思える。これは⃗
yをより簡素に説明することで、そ22∥⃗x∥2 ≡ |x1 | + · · · + |xn |れを正則化という。
• 0 次ノルムは注意が必要である。0 次ノルムは ⃗x の非ゼ
尚、0 次ノルム最小化でスパースな ⃗x を完全に再構成
ロ要素の数を意味し、次式で書かれる。
できることは証明[12,13]
されている。
n∑∥x∥0 =|xi |0 = |x1 |0 + · · · + |xn |0
しかし、0 次ノルム最小化を解くのに必要な計算量は、組
i=1
み合わせ爆発が起こるため、ω
⃗ の個数 N に応じて指数的にしかし、非ゼロ値 a (̸= 0) では、a0 = 1 であるのに対し、 増大し、高次元の問題に応用するのは現実的ではない。
本来数学的には 00 は不定である。厳密にはそうである
が、00 = 0 として定義されたものが 0 次ノルムである。 3.3.2 L1 正則化: LASSO 法
Tibshirani[14] は線形回帰の問題に 1 次ノルムの項を加つまり 0 次ノルムは、⃗
x の非ゼロ要素数を意味する。
える “LASSO 法” (Least Absolute Shrinkage and Selection
3.1.2 arg max (min)Operator) という手法を提案した。
arg max (arg min) は、引数 x を必要とする関数 f (x) が最L1 正則化: LASSO 法
大値 (最小値) となる時の{ x を返す関数で、}LASSO 法で取り扱う問題は(44) 式と書ける。
(
)
−
→ 2
arg max f (x) ≡ x ∈ S : f (x) = max f (y)
b
ω
⃗
=
arg
min
∥⃗
y
−
X
ω
⃗
∥
+
α∥⃗
ω
∥
(44)
1
2
y∈S
x∈S
ω
⃗
{
}
• α = 0 : 最小二乗法と同じ解を与えるが、それはノイ
arg min f (x) ≡ x ∈ S : f (x) = min f (y)y∈S
ズまで合わせに行く過適合 (over-fitting) な解である。
x∈S
で定義される。
• 0 < α < ∞ : 適切な規範でスパース解を得る。
⃗x の横軸に対し ⃗y の測定データ・セットが得られた場合、 • α = ∞ : ∥⃗ω ∥1 = 0 が要請される。つまり非ゼロ要素⃗ とする
そのデータ・セットを説明する物理モデルを f (⃗x; θ)はなく、データは全てノイズと解釈される。
⃗ の平均二乗と、最小二乗法は測定値 ⃗y と物理モデル f (⃗x; θ)《講義ノート》→−
LASSO 法は、 X の列ベクトルの中で不要なものは係数ω⃗ をゼロにすることで過適合 (over-fitting) を防ぎ、モデルの予測誤差を小さくすることを目的としている。また、ω
⃗の中でデータを説明するのに必要な変数だけを選択する手法 (変数選択) としても使われる。
1 次ノルム最小化の解 = 0 次ノルム最小化の解
さらに、2000 年代半ばには Candes & Tao[15] やDonoho & Tanner[16] が、M < N の問題 (不良設定問題) で、ある条件下においては 1 次ノルム最小化に
よって得られる解が 0 次ノルム最小化で得られる解と同じになることを数学的に証明した。
1 次ノルム最小化は従来からよく知られた最適化アルゴ
リズムである「線形計画法」で解くことができ、計算量の面では 0 次ノルム最小化より有利になる。

物性研究・電子版 Vol. 8, No. 1, 081201（2020年2月号）• α∥⃗ω ∥1 :
α ≥ 0 であるから、原点を極小値として上に開いた V字型関数である。この「原点を極小値」にとる特徴がα → ∞ で ω1 → 0 を導き、パラメータ ω1 をゼロへ誘導しスパース化させる。
図 19-(a)∼(e) の赤破線で示した様に、(−→ 2)−1.4 = arg min ∥⃗y − X ω⃗ ∥2
ω1

で極小値 (図中の min1 ) を取るのに対し、青破線で示した
−
→ 2
様に ∥⃗
y−Xω
⃗ ∥2 + α∥⃗ω ∥1 がとる極小値 x̂(図中の min2 ) は
α に依存する。これらのの α 依存性を図 19-(f) に示した。
図 19-(f) より、α > 15 では ωˆ1 ≈ 0 となり、唯一非ゼロであった ω1 もゼロに漸近する。その結果、ω
⃗ = 0 となり、測定で得られた ⃗y = (−2, 3) が残差として残ることになり、主成分を発見できずにデータは全てノイズである解析となる。

3.3.3 L1 正則化で何故スパースな解が得られるか
直感的に理解するために、観測量 ⃗y が 1 つ (1 次元: M = 1) 3.4 時系列データの周期解析への応用
この節では図 20 の青線で示したようなでパラメータ ω⃗ が 2 次元 (N = 2) の場合について考える。
N∑
この場合 (42) 式の連立方程式は次式となる。
 
y(t)
=
(aj cos 2πfj t + bj sin 2πfj t) + n(t)
(
) ω
1
j=1
(y1 ) = X11 X12   → y1 = X11 ω1 + X12 ω2 ,でモデル化出来るような時系列データを対象とする。
ω2
図 20 のデータは、f1 = 0.8 Hz, f2 = 2.4 Hz, bf1 = 1.0,
もちろん方程式が 1 つで、未定係数が 2 つ (ω1 , ω2 ) なので、
af
2 = 0.8 で、σnoise = 0.2 の標準偏差を持つ正規分布ノイ
ユニークな解は得られない。
ズ N (0, σnoise ) で、
1 次ノルム最小化の模式図
y(t) = bf1 sin 2πf1 t + af2 cos 2πf2 t + N (0, σnoise )方程式は、で合成したものである。
ω2
X11
y1
ω2 =
ω1 +
X12
X12
3.4.1 フーリエ変換基底の用に直線 (図 18 の赤線) を表
このような振動性の時系列データから振動周波数を推定し、その直線上の任意の点でする場合、フーリエ変換がよく用いられる。しかしフーリ−→⃗y = X ω⃗ の関係を満たす。
エ変換は完全性が担保される結果、データに必ず重畳するω1
一方図 18 の青線は内側からノイズまでも周波数分解してしまう。
順に 1 次ノルムが大きい場合
ここではフーリエ変換の基底関数を用いて SpM を試みを示しているが、L1 正則化はる。時系列データが ∆t の等時間ステップで M 個のデータ
赤い直線の内、1 次ノルムが最 図 18: 赤い直線が方程式を点として計測されたとすると、フーリエ変換の周波数分解
小のものを得る手法である。こ 満たす線で、青い線は 1 次
能 ∆f は、∆f ≡ 1/(M ∆t) である。ここで fn ≡ n∆f とすの例では一番内側のひし形がそ ノルム一定の線である。
ると、N = M/2 としてフーリエ変換は
N
れである。
∑
図から直線と ω1 軸と ω2 軸と交わる点は 2 箇所あるこ(an cos 2πfn t + bn sin 2πfn t)(45)y(t) =とが変わる。それらの交点ではともに |⃗ω |0 = 1 である。
n=1
L1 正則化は、その |ω|0 = 1 のスパース解のうち、|⃗
ω |1 がである。y(t) の時系列データを ⃗y ≡ {y1 , y2 , · · · , yM } とし、より小さい ω2 軸との交点が選ばれる。
ω⃗ ≡ {a1 , a2 , · · · , aN , b1 , b2 , · · · , bN } とすると、(45) 式は、
 cos 2πf t · · · cos 2πf t sin 2πf t · · · sin 2πf t 
この結果は、主要パラメータ抽出の観点から考えると、
1 1
1 1
N 1
N 1
 cos 2πf1 t2 · · · cos 2πfN t2 sin 2πf1 t2 · · · sin 2πfN t2 
ω1 よりも ω2 のほうが、⃗y を再現するには重要であるこ


..
..
..
..


..
..
とを意味している。


.
.
.
.
.
.


−
→ 

X
≡
.
.
.
.
さらに (44) 式について理解を深めておく。不良設定問題
..
..


.
.
.
.


.
.
.
.
.
.


でなくなるが、
 
 


.
.
.
.
.
.


..
..
..
..
..
..
(
)
−2
2
−
→
⃗y =   , X =   , ω⃗ = ω1cos 2πf1 tM · · · cos 2πfN tM sin 2πf1 tM · · · sin 2πfN tM
−
→
−
→
3
−1
を用いて、⃗
y = Xω
⃗ と書ける。ここで X を図 21 に示しとしよう。この場合、た。 この縦一列が各周波数の基底関数である。

−2 − 2ω1−→
3.4.2 α に依存した非ゼロ要素の変化
,⃗y − X ω⃗ →
3 + ω1
図 22 に ω⃗ 係数 (= {a1 , a2 , · · · , aN , b1 , b2 , · · · , bN }) の非ゼロ要素がα の減少とともに分岐する様子を示した。
−→ 2∥⃗y − X ω
⃗ ∥2 = (−2 − 2ω1 )2 + (3 + ω1 )2
α の大きい右側では適切に 2 つの振動モードが抽出されているのに対し、α の小さい左側では、小さな非ゼロ要素= 5ω12 + 14ω1 + 13が多数分岐し始めることが分かる。その領域は、データに青z}|{重畳したノイズ成分をフーリエ基底で再現しようとする過−→ 2∥⃗y − X ω
⃗ ∥2 + α∥⃗ω ∥1 = 5ω12 + 14ω1 + 13 + α|ω1 |
適合領域である。
{z} | {z }|図 23 は、α を横軸 (対数) にして、○は ω⃗ の非ゼロ要素数
赤
緑
∥⃗
ω
∥
の変化を、●は得られたスパース解ω⃗ˆ による再現デーである。これらをグラフに示したものが図 19 である。
0−→ˆ
ここでこれらのグラフ特徴をまとめておく。
タX ω⃗ と解析対象データ ⃗y との残差二乗平均平方根 (RMSD:−→ 2⃗ ∥2 :• ∥⃗y − X ω
Root Mean-Square of Deviations) の変化を示した。水平破線
ω1 に関して上に開いた放物線関数であり、その極小値は はデータに重畳させた正規ノイズの標準偏差である。
必ず 0 以上である。

《講義ノート》

物性研究・電子版 Vol. 8, No. 1, 081201（2020年2月号）
α=1.,
min1 =-1.4,
min2 =-1.320
5ω21 +14ω1 +13

15

Out[32]=α|ω1 |

- 2.0

- 1.5

- 1.0

- 0.5

α=3.,
min1 =-1.4,
min2 =-1.1

0.0

ω1

0.5

(b)

- 2.0

- 1.5

- 1.0

- 0.5

20
5ω21 +14ω1 +13

α|ω1 |

- 1.0

- 0.5

5ω21 +14ω1 +13+α|ω1 |

α|ω1 |
5ω21 +14ω1 +13+α|ω1 |

10
5

0.0

ω1

0.5

(d)

- 2.0

- 1.5

α=15.,
min1 =-1.4,
min2 =-1.42706 × 10-8

- 1.0

- 0.5

ω1

0.0

0.5

%"%

()

%"$

,,-( (. w,,$(/ a,,w,,!

!#w (a*0(
1

$

!$

%"#
*a(1w

20

5ω21 +14ω1 +13)(

15

Out[48]=5ω21 +14ω1 +13

15

Out[44]=α|ω1 |

&

%"&

'

!"%

5ω21 +14ω1 +13+α|ω1 |

10

!%

%"'

$

- 1.5

0.5

20

5

- 2.0

ω1

0.0

α=10.,
min1 =-1.4,
min2 =-0.4

10

(c)

5ω21 +14ω1 +13+α|ω1 |5

15

Out[40]=α|ω1 |

105(a)5ω21 +14ω1 +13

15

Out[36]=5ω21 +14ω1 +13+α|ω1 |

10

20

#

- 2.0

(e)

- 1.5

- 1.0

- 0.5

!

!"$

5

,,w,,a /($,,w .( (-,,

α=0.3,
min1 =-1.4,
min2 =-1.37

$

!"#

0.0

ω1

0.5

!%

(f)

#

!%

+

!%

$

!%

!

%

!%

%

!

$

!%

!%

a

図 19: (a)∼(e) α = 0.3 ∼ 15。(f) α 依存性

92

66

−7

−6

−5

L0 norm
28
2

2

2

−2

−1

1.0

97

0.0

−2

0.2

−1

Y
0

Coefficients
Coefficients
0.4
0.6

0.8

1

2

Y
Ym
noise

0

2

4

6

8

10

time

図 22: ω⃗ の α に依存した分岐の変化。横軸は α を対数目盛りで示してある。

50

10

図 20: 時系列データの周期解析対象のデータ−4−3ln(lambda)

RMSD

1.0

40

8

0.8

Znorm

30
Znorm
0

−1.0

0.2

2

10

−0.5

0.4

4

20

0.0

RMSD
0.6

time

6

0.5

20

40
60
80
cos(2*pi*f*t), sin(2*pi*f*t)100

図 21: 線形写像行列のイメージプロット

−3.0

−2.5

−2.0
−1.5
−1.0
log10(lambda)−0.5

図 23: 0 次ノルムと RMSD の α による変化のグラフ

《講義ノート》

物性研究・電子版 Vol. 8, No. 1, 081201（2020年2月号）

log10 (α) > −1.5 では 2 つの振動モードが抽出され、α の減少とともに RMSD が減少して再現性が向上していることがわかる。一方、log10 (α) < −1.5 ではノイズを再現するた
め 0 次ノルムが急激に増加し、RMSD が重畳させた正規ノ
イズの標準偏差を下回る過適合になっていることが分かる。
3.4.3 交差検定
SpM では、データを説明する成分の冗長性とデータの再現性両方を考慮した適切な情報量規準にしたがって α を決定する。α を決める方法として、交差検定 [17] などもよく用いられる。

3.5

EXAFS 解析への応用

Photoelectron wave
XrayPhotoelectron wave
Xra
y

hu

E

E

EXAFS

EXAFS振動
EXAFS振動

hu

EXAFS
E0

md

y

ray

ra
X-

X-

md

E0

図 24: 広域 X 線吸収微細構造 (EXAFS)LASSO 法などの L1 正則化による SpM は、抽出したい
特徴量 ω⃗ が少数 (スパース) である場合、従来解析法を凌駕しうる強力な解析法である。
赤井等はその優位性に注目し、放射光を用いて計測される
広域 X 線吸収微細構造 (EXAFS: Extended X-ray Absorption
Fine Structure) の解析に L1 正則化を適用 [18] した。
図 24 に示した様に、EXAFS は、X 線吸収端以上への電子励起による光電子波の近接原子による散乱と干渉によって、吸収スペクトルに現れる振動構造である。EXAFS では、X 線吸収端エネルギーで選択された注目原子近傍の原子による散乱が重要であるが、本来近接原子の配位は、物質の化学結合や構造を反映して、距離に対してスパースである。赤井等はそのスパース性に注目した。

0.12

CVE(λ,σDW)
0.155

0.145
^ DW = 8.24 × 10−2
σ

0.140
0.135
^
λ = 3.585 × 10−3

0.04

°)
Debye Waller σDW (A
0.06
0.08
0.10

0.150

10−5

10−4
10−3
λ (logarithmic scale)

0.130
0.125
10−2

図 25: 交差検定誤差のヒートマップ

EXAFS では、注目原子近傍の近距離構造の動径分布関
数 (RDF: Radial Distribution Function) を解明する目的で解析されるが、近接配位原子の構造ゆらぎや可動性を表すデバイ・ワラー因子を推定するには、本来未解明の近距離構造を仮定する必要があった。それに対し、L1 正則化を用いた解析では、一切の事前仮定を必要とせず、データだけから RDF とデバイ・ワーラー因子の推定に成功した。その例を図 25 に示した。この横軸はスパース性を制御する α を対数で示し、縦軸はデバイ・ワラー因子で、交差検定誤差をヒートマップで示したものである。交差検定誤差の最小化
によってデバイ・ワーラー因子が推定可能である。

謝辞
本文の一部は、JST, CREST, JPMJCR1861 の支援を受けたものである。

参考文献[1] T. Bayes and R. Price. An Essay towards Solving a Problem in the Doctrine of Chances. Phil. Trans. Roy. Soc.,
53, pp. 370–418, (1763).
[2] K. Nagata, S. Sugita, and M. Okada. Bayesian spectral
deconvolution with the exchange Monte Carlo method.
Neural Networks, 28, p. 82, (2012).
[3] P. J. Green. Reversible jump Markov chain montecarlo computation and Bayesian model determination.
Biometrika, 82, p. 711, (1995).
[4] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth,A. H. Teller, and E. Teller. Equation of state calculations by fast computing machines. J. Chem. Phys., 21,
pp. 1087–1092, (1953).
[5] C. Andrieu, N. D. Freitas, A. Doucet, and M. I. Jordan. An
Introduction to MCMC for Machine Learning. Machine
Learning, 50, pp. 5–43, (2003).
[6] S. Aihara, M. Hamamoto, K. Iwamitsu, M. Okada, andI. Akai. High precision modeling of a damped oscillationin coherent phonon signals by Bayesian inference. AIP
Adv., 7, p. 045107, (2017).
[7] R. H. Swendsen and J-S. Wang. Replica Monte Carlo
Simulation of Spin-Glasses. Phys. Rev. Lett., 57, pp.
2607–2609, (1986).
[8] K. Hukushima and K. Nemoto. Exchange Monte Carlo
Method and Application to Spin Glass Simulations. J.
Phys. Soc. Jpn., 65, pp. 1604–1608, (1996).
[9] Y. Mototake, M. Mizumaki, I. Akai, and M. Okada.
Bayesian Hamiltonian Selection in X-ray Photoelectron
Spectroscopy. J. Phys. Soc. Jpn., 88, p. 034004, (2019).
[10] 岡田真人. 新学術領域研究「スパースモデリングの深化と高次元データ駆動科学の創成」, 2014. http://sparse-modeling.jp/.
[11] 植 村 誠.
天 文 屋 の た め の HOW TO ス パ ー ス
モ デ リ ン グ, 連 立 1 次 方 程 式 を 解 く, 2016.
http://home.hiroshima-u.ac.jp/
uemuram/How2SM/01_demo.pdf.
[12] 和田山正. 圧縮センシングにおける完全再現十分条件について. 日本神経回路学会誌, 17, pp. 63–69, (2010).
[13] 田中利幸. 圧縮センシングの数理. IEICE Fundamentals Review, 4, pp. 39–47, (2010).
[14] R. Tibshirani. Regression shrinkage and selection via thelasso. J. Roy. Stat. Soc. B, 58, pp. 267–288, (1996).
[15] E. J. Candes and T. Tao. Near-optimal signal recovery
from random projections: Universal encoding strategies?
IEEE Transactions, Information Theory, 52, pp. 5406–
5425, (2006).
[16] D. L. Donoho. Compressed sensing. IEEE Transactions,
Information Theory, 52, pp. 1289–1306, (2006).
[17] K. P. Murphy. Machine Learning: A Probabilistic Perspective (Adaptive Computation and Machine Learning
series). The MIT press, 2012. ISBN 0262018020.
[18] I. Akai, K. Iwamitsu, Y. Igarashi, M. Okada, H. Setoyama, T. Okajima, and Y. Hirai. Sparse Modeling of
an Extended X-Ray Absorption Fine-Structure Spectrum
Based on a Single-Scattering Formalism. J. Phys. Soc.
Jpn., 87, p. 074003, (2018).