パラメータ推定解説

数理モデルに対するパラメータ推定長山 雅晴
北海道大学 電子科学研究所
附属社会創造数学研究センター 人間数理研究分野

松永哲
北海道大学 理学院 数学専攻

自己紹介
所属：北海道大学電子科学研究所
附属社会創造数学研究センター 人間数理研究分野
専門：応用数学
キーワード：
現象を数式であらわす数理モデリングとその数値シミュレーション
化学反応や自然界に現れるパターン形成の数理解析
数理モデルの解析をおこなうための数値分岐計算
1988年：徳島大学・総合科学部，卒業
1992年：広島大学・理学研究科数学専攻，博士前期課程，修了
1994年：東京大学・数理科学研究科，博士課程後期，修了
1998年：龍谷大学・ハイテクリサーチセンター，博士研究員
1999年：京都大学・数理解析研究所，助手
2004年：金沢大学・大学院自然科学研究科，助教授
2009年：金沢大学・理工研究域数物科学域，教授
2012年：現在

数理モデリング
自然科学・人間社会における諸問題の解決
自
走
粒
子
運
動

反
応
拡
散
系

理
学

燃
焼
合
成

焼
結
過
程

工
学

細
胞
核
運
動

分
化
の
波

慢
性
白
血
病

生
命
科
学

ガ
形ン
成腫
問瘍
題の

医
学

皮
膚
応疾
用患
へ
の

渋
滞
現
象

社
科会
学・
環
境

数理モデリング（諸分野と数学を繋ぐ接着剤）
数学・数理科学

汚
水
処
理
問
題

講義の予定1. 数理モデリングとは何か？
2. 数理モデリングの評価
1.
2.

3.

パラメータ推定（1点推定）基礎編
1.
2.

4.

定性的な評価
定量的な評価解説
演習

ベイズ推定を用いたパラメータ推定
1.
2.

解説
演習

数理モデリングとは何か？

自然現象（生命現象）
理学
３つの手法

自然現象を科学的に理解する
実験科学

理論科学

計算科学

数理モデリング
数理モデリングは理論的な理解を与えるための道具である．
現象の定性的な再現は必要条件である．
理論的方法によって未知の部分が理解できなければ意味がない．
数理モデルに正解はない（40点〜80点程度）
理解できれば，実験に対する予測ができる．
よい数理モデルができれば，数学の問題へと定式化ができる．

数理モデルとは何か？
数理モデルとは現象の数式による模倣と抽象化である．
抽象化し現象を模倣することによって，「現象の法則」を明らかにする道具

モデル化

BZ反応（現象）

抽象化

モデルから得られた結果
法則の発見抽象化したモデル

どのようなモデルを考えればよいか？（立場によるが）現象の定性的性質を保持してる程度まで抽象化した数理モデルを考えたい．
そのような数理モデルは現象の法則を記述してハズ（信念！）数理モデルは現象を数学で取り扱えるようにする道具ともいえる．

数理モデルはどのように作るか？

l og

cn (an )cn+1)
na) a f(
数理モデル（漸化式，微分方程式）nnbs( b nc n +1⇢oc nan+1 an = f (an , bn ),h(a a an cn g(an , bn )n, bbn+1 bn = g(an , bn ).
bn+1n)
数理モデルは)bn)n どのように作るの？
n
=
0,
1,
2,
·
·
·
a
n
(
b
s
(
o
g
c
o
l

sin

an
b
n

)
(a n

私とっての数理モデルを作る基本原理
重要

単位時間当たりの（量，もの）の変化量＝ 流入量 ー 流出量＝ 入ってくる（もの）ー 出て行く（もの）＝ 生成されるもの ー 消費されるもの(もの) ＝ 化学物質濃度，密度，温度，力

etc.

数理モデルはこの基本原理によって記述されていると思ってよい．
「流入するもの」と「流出するもの」をいかにして見つけるのかが重要！
いいかえると「“ルール”をいかにして見つけ出すか」がポイント

数理モデリングに必要な道具

力学系の基礎知識
常微分方程式系の平衡解の線形化安定性解析
相図
簡単な分岐現象
数値計算法の基礎（水藤先生，斎藤先生）
常微分方程式の数値計算法
熱方程式の数値計算法（できれば）反応拡散系の数値計算
現象に応じて，物理学，化学，生物学，その他

数理モデルの構成

・知りたい現象（実験結果）
・自然科学言語を用いて現象の説明をする
・数式で表現（微分方程式，漸化式など）
・無次元化［変数，定数の持つ単位を「1」にする］

No・数値計算による現象再現の評価（数理モデルの評価は？）
Yes

・数理モデルを用いた実験の予測
・実験系の提案

現象再現しただけでは数理モデルを評価した
ことにならない．

数理モデルの評価とは現象の数理モデリングを評価する方法は？
現象を再現しているだけでは不十分（似て非なる数理モデルを作っているだけの可能性あり）
現象を矛盾なく説明できているならば新しい現象を予測できるのではないか？
モデルからの予測と実験によるその再現性から数理モデルを評価できないか？
定量性の観点から数理モデルを評価することはできないか？
定量性はどのようにあたえることができるのか？
定量的再現性の視点から数理モデルを評価してみたい…

似て非なるもの（リングパターン）
生物実験

数値実験

Wakita et al., 1994.

リングパターンの成長様式が異なる

Ismael Ráfols

Mimura et al., 2000.

似て非なるもの（スパイラルパターン）
単安定反応拡散系モデル

BZ反応モデル

燃焼反応モデル

双安定反応拡散系モデル

BZ反応モデル

✓
8
@u
1
a u
>
>
< @t = Du u + " f v a + u + u(1@v>>= Dv v + u v.
:
@t

単安定興奮系◆u) ,発熱反応拡散系
供給

k(T)A + O → P + Heat
対流による放熱15発熱反応モデルのダイナミクス

単安定興奮系

双安定系反応拡散方程式

似て非なるもの（スパイラルパターン）
単安定反応拡散系モデル

双安定反応拡散系モデル

BZ反応モデル

燃焼反応モデルパターンが定性的に似ているからといってメカニズムが同じであるという保証は全くない

数理モデルの評価とは現象の数理モデリングを評価する方法は？
現象を再現しているだけでは不十分（似て非なる数理モデルを作っているだけの可能性あり）
現象を矛盾なく説明できているならば新しい現象を予測できるのではないか？
モデルからの予測と実験によるその再現性から数理モデルを評価できないか？
定量性の観点から数理モデルを評価することはできないか？
定量性はどのようにあたえることができるのか？
定量的再現性の視点から数理モデルを評価してみたい…

定量的再現性のある数理モデリングを目指して
すでに実験データから偏微分方程式モデルを構築する研究が行われている

現象
実験

機械学習

可能なのか？
可能な場合の十分条件が
あるはず．

未来予測
制御

定量性を確保できれば検証

数理モデリング
数理モデルのパラメータを決定できるための意味のある実験データ（データ構造）は何か？
明らかにすることが必要

定量性のある数理モデリング
実験データが与えられていると仮定する．
このときのある数理モデリングを行うためには2つの問題がある（1）数理モデルは現象を記述しているのか？
（2）数理モデルはデータを再現しているが，定量性はあるのか？
(1) に関しては現象の数理モデリングの定性的評価を行う必要がある．
ここでは (2) について評価方法を考えてみる．
(2) についての問題点：与えられたデータから(2.1) モデル方程式のパラメータを見つけることができるのか？
(2.2) モデル方程式のパラメータを決定できるのか？
←どのようなデータならばパラメータを決定できるのか？

(2.2) は谷口（やぐち）さん（神戸大学）のCRESTで研究されています．

問題設定（パラメータ推定）
数理モデリングは職人芸なので対象としない．
（本当は数理モデリングの話もしたいが…）
定性的に評価された数理モデルを仮定する．
すなわち，モデル方程式<latexit sha1_base64="YzSs3JDpJQJT1XvoSkoLZO+iAGk=">AAADU3ichVNNaxNBGH6T+FFXa2O9iF4Wk0oKkk7qoaII1SJ4KbZp0wYyYZndTNKhsx/MToLrsn/Ag1cPnhQ8iD/Di3+gSH+CeGzBgx58d7JQbFFnmZlnnnmfd56Zl3UjKWJNyGGpXDl3/sLFmUvW5SuzV+eq1+Z34nCsPN7xQhmqrstiLkXAO1poybuR4sx3Jd9199fy/d0JV7EIg22dRLzvs1EghsJjGimn+rpet6nkQ01Ti7p8JIKUKcWSLJWZRQcijiRLYp1IbtOhYl46sKnrjzOcdWY/yhfDRs48zCFbvEtp757v909rMaKhF6eCsUOaFuXBoDjKokqM9nTTqtedao00iWn2WdAqQA2KthFWvwKFAYTgwRh84BCARiyBQYxfD1pAIEKuDylyCpEw+xwysFA7xiiOEQzZfRxHuOoVbIDrPGds1B6eIrErVNqwQA7IR3JEvpBP5Bv59ddcqcmRe0lwdqdaHjlzr25s/fivysdZw96J6p+eNQzhvvEq0HtkmPwW3lQ/efnmaOtBeyG9Q96T7+j/HTkkn/EGweTY+7DJ229N9hNtCkvQMWcpHJfMe4ywJ8YdQ2YbnkIXHFiHx7AGbXiO3Hrh2THIQ20ITePkBWRY4Nbpcp4FO8vNFuLN5drqk6LUM3ALbkMD67kCq/AMNtCZB8elm6VaqV4+KP+s4F8yDS2XCs11+KNVZn8DZxvOBQ==</latexit>

に対して， 初期値を込めたパラメータセットを与えて数値計算から数値データを求め，この数値データ（にノイズを加えたデータ）を与えられたデータとする．

問題
このとき，数値データを近似するパラメータセット

モデル方程式から求めることができるか？（推定できるか？）を

最適化問題として定式化
モデル方程式<latexit sha1_base64="YzSs3JDpJQJT1XvoSkoLZO+iAGk=">AAADU3ichVNNaxNBGH6T+FFXa2O9iF4Wk0oKkk7qoaII1SJ4KbZp0wYyYZndTNKhsx/MToLrsn/Ag1cPnhQ8iD/Di3+gSH+CeGzBgx58d7JQbFFnmZlnnnmfd56Zl3UjKWJNyGGpXDl3/sLFmUvW5SuzV+eq1+Z34nCsPN7xQhmqrstiLkXAO1poybuR4sx3Jd9199fy/d0JV7EIg22dRLzvs1EghsJjGimn+rpet6nkQ01Ti7p8JIKUKcWSLJWZRQcijiRLYp1IbtOhYl46sKnrjzOcdWY/yhfDRs48zCFbvEtp757v909rMaKhF6eCsUOaFuXBoDjKokqM9nTTqtedao00iWn2WdAqQA2KthFWvwKFAYTgwRh84BCARiyBQYxfD1pAIEKuDylyCpEw+xwysFA7xiiOEQzZfRxHuOoVbIDrPGds1B6eIrErVNqwQA7IR3JEvpBP5Bv59ddcqcmRe0lwdqdaHjlzr25s/fivysdZw96J6p+eNQzhvvEq0HtkmPwW3lQ/efnmaOtBeyG9Q96T7+j/HTkkn/EGweTY+7DJ229N9hNtCkvQMWcpHJfMe4ywJ8YdQ2YbnkIXHFiHx7AGbXiO3Hrh2THIQ20ITePkBWRY4Nbpcp4FO8vNFuLN5drqk6LUM3ALbkMD67kCq/AMNtCZB8elm6VaqV4+KP+s4F8yDS2XCs11+KNVZn8DZxvOBQ==</latexit>

に対して
パラメータ

を与えて求まる数値解をとする

このとき，次の目的関数を与える：

この目的関数の最小化から正しくパラメータを推定することができるか？

パラメータ推定方法

をパラメータとして，次式を最小化する

勾配法(1)

を求める

目的関数

数値計算結果と実験値の差の2乗ノルム

観測誤差の過剰学習(オーバーフィッティング)を防ぐ役割(: 正則化パラメータ)

パラメータ推定方法

適当なパラメータの初期点下げるように点列パラメータ

勾配法(2)

から出発して目的関数の値をを生成し，最適解に収束させる(点推定).
から次のパラメータに以下のように更新する

ステップ幅
探索方向として勾配

探索方向を用いる方法を最急降下法と呼ぶ．

局所最適解に陥りやすい（欠点）最急降下法1.

を初期値として与える. k = 0とする

2．次の三つの条件を同時に満たした場合に停止する

3．探索方向を dk =<latexit sha1_base64="HQntFH8IutQIuoRExRb2WDej8qo=">AAAC9nichVHPaxNBFP66aq1R26gXwctiUqkHk0lFFEGoFsFLsU2bNtAty+xkki6Z/cHuJDQu/Qe8iwdBUCgiXrx7VMF/QLB/gnis4MWDbycLRYv6hpn35pv3vflmnhcrP9WM7U9Yx46fmDw5dap0+szZ6ZnyufPraTRIhGyJSEVJ2+OpVH4oW9rXSrbjRPLAU3LD6y/m5xtDmaR+FK7pUSy3At4L/a4vuCbILder1ZLjBR23b9+xrzkh9xS3u3OOF6lOOgrIZY7elprvuv2rpWrVLVdYjRmzjwaNIqigsOWo/AUOOoggMEAAiRCaYgWOlMYmGmCICdtCRlhCkW/OJXZRIu6AsiRlcEL7tPZot1mgIe3zmqlhC7pF0UyIaWOWfWav2QH7xN6wr+znX2tlpkauZUTeG3Nl7M48vrj647+sgLzG9iHrn5o1urhltPqkPTZI/gox5g8fPT1Yvd2cza6wl+wb6X/B9tl7ekE4/C72VmTzmal+yM1QR8vcldBaN//Rozky6jgha7iPNlws4S4W0cRDwpYKza6JBHEj1IySnfzPqcWNPxt6NFifrzWu126szFcW7hXNnsIlXMYcdfQmFvAAy6RN4Ane4QM+WjvWc2vPejVOtSYKzgX8ZtbbX3jTrfQ=</latexit>

rf (✓ k ) と定める(↵) = f (✓ k + ↵dk ) に対する直線探索により ↵k を決定

4.

<latexit sha1_base64="CH4tHDYUhajD9Ax54KE6ZmiNt/U=">AAAC03ichVFBSxtBGH1uW2ujNbFeBC/SJKWnZKIUS0/RUPASqonRQJIus+skWbLZXXY3wbh4kZ4KPbfQk0IPpT+jl/4BwfwE6dFCLz347WQhaNDOMDNv3nzvmzfzaY5peD5jwynlwcNH049nnsRm557OxxMLz/Y8u+fqoqLbpu1WNe4J07BExTd8U1QdV/CuZop9rVMIz/f7wvUM29r1B45odHnLMpqGzn2i3qdSsTo3nTZXO7FUSk0kWYbJtjIJchFIImrbduICdRzAho4euhCw4BM2weFRryEHBoe4BgLiXEKGPBc4Roy0PYoSFMGJ7dDcol0tYi3ahzk9qdbpFpOGS8oVpNk5+86u2C/2g12yf3fmCmSO0MuAVm2kFY4a/7hU/vtfVZdWH+2x6l7PPpp4Lb0a5N2RTPgKfaTvH32+Kr8ppYMX7Iz9Jv+nbMh+0gus/h/9244ofZXZx9oAWVTkXS7NWfkfLRoD6Y4Ts4u3qEJFERsooIR3xBUjz6pEOmltZKSTw/DPqcS52wWdBHurmdxa5tXOajK/GRV7Bst4jpdU0XXksYVt8hZm/4JTnCkVJVBOlA+jUGUq0iziRlM+XQOW2aCq</latexit><latexit sha1_base64="9LSK713fpshL1nzyKuBBrAe4kg0=">AAADJHichVHPaxNBFP66/miNPxr1InhZDJVUIZ1UiiII1SJ4KbZp0wa6YZndTJIh+6u7k2Bc8g/4D3jwpKAg/hm9ePAmBXvxXrwZwYsH304Wqhb1DTPvzTfve/PNPCfyZKIYO5gyTpw8dXp65kzh7LnzF2aLFy9tJWE/dkXdDb0wbjg8EZ4MRF1J5YlGFAvuO57Ydnor2fn2QMSJDINNNYxE0+edQLalyxVBdrFpWo7oyCAVu30N3RiZBdOKutIsW9yLunz+XrtsOaHXSoY+udRSXaH4yO6ZN81JBpXwW3ZvnngiaP1ayS6WWIVpM48H1TwoIbe1sPgJFloI4aIPHwIBFMUeOBIaO6iCISKsiZSwmCKpzwVGKBC3T1mCMjihPVo7tNvJ0YD2Wc1Es126xaMZE9PEHNtnb9mYvWfv2CH78ddaqa6RaRmSdyZcEdmzz65sfP8vyyev0D1i/VOzQht3tFZJ2iONZK9wJ/zB0+fjjbu1ufQ6e8W+kP6X7IDt0QuCwTf39bqovdDVj7gpFlDXd8W0Luj/6NAcanWckE08RAM2VnEfK6jhMWGruWZbRy5xQ1S0kifZn1OLq3829HiwtVip3qosrS+Wlh/kzZ7BVVxDmTp6G8t4hDXS5mIPh/iKsfHG+GB8NPYnqcZUzrmM38z4/BPJcsFK</latexit>

5.

✓ k+1 = ✓ k + ↵k dk と更新<latexit sha1_base64="NkYYnEhhqyRRrDiPoixL+S4I+aE=">AAADEnichVG/axRBFP6y/ornj5ymEWwW7yJC4DIbEUUIRIMgSDC55JKDXFhm9yZ3y87+YHfu8FzuH7AWLKwULMTK1k7SWNiGmE5LsYxgY+HbuYWg0fiGmXnzvfe9+WaeE0svVYztjRnHjp84eWr8dOnM2XPnJ8oXLq6lUS9xRcONZJQ0HZ4K6YWioTwlRTNOBA8cKdYdfyGPr/dFknpRuKoGsdgMeCf0tjyXK4Ls8oNqtdRyItlOBwFtWUt1heJDO/OnraE5Z/49NjSnzRaXcZfbPqUEbdsvVat2ucJqTJt52LEKp4LClqLyJ7TQRgQXPQQQCKHIl+BIaWzAAkNM2CYywhLyPB0XGKJE3B5lCcrghPq0dui0UaAhnfOaqWa7dIukmRDTxBTbYa/ZPvvA3rCv7Oc/a2W6Rq5lQLsz4orYnnhyaeXHf1kB7QrdA9aRmhW2cEtr9Uh7rJH8Fe6I33/8bH/ldn0qu8pesm+k/wXbY9v0grD/3X21LOrPdfUDboYZNPRdCa0z+j86NAdaHSdkFffQhI1F3MEC6nhI2GKh2daeS9wINa3kUf7n1GLrz4YedtZma9b12o3l2cr83aLZ47iMK7hGHb2JedzHEmlz8Ra7+IwvxlPjnfHe2B6lGmMFZxK/mfHxF4X5ung=</latexit>

6. k = k+1としてステップ2に戻る
金森敬文, “機械学習のための連続最適化” (講談社, 2016)

αの選択方法
方針1 まじめに1次元の最適化問題を解く→を固定するとはαに関する
1次元の最適化問題になるのでこれを真面目に解く
・二分法
・黄金分割探索方針2 大体でいいや→減少すればいいので，大枠で評価してしまう・Wolf の基準
・etc.

参考文献：
金森敬文，鈴木大慈，竹内一郎，佐藤一誠著，
機械学習のための連続最適化，講談社

直線探索↵ に関する1変数関数(↵) = f (✓ k + ↵dk ) (↵ > 0)<latexit sha1_base64="lR1tPdba67cuLJznjCrGeU9dE7s=">AAAC0XichVFNSxtBGH5cP2rjR6K9CL2IieIpTiyiePKDQi+iJiYGjITZcZIsbnaX3U0whoJ4bO/twYOt4EH8GV76Bwr6E8RjBC8efHeyICrqDDPzzDPv884z8+qOaXg+Y1cdWmdXd8+H3o+Rvv6BwWhsaDjn2TVXyKywTdvN69yTpmHJrG/4psw7ruRV3ZSb+u5ycL5Zl65n2NaG33DkdpWXLaNkCO4TVUgkIgVuOhUeSSSKsThLMtVGX4JUCOII25odu0QBO7AhUEMVEhZ8wiY4POpbSIHBIW4bTeJcQoY6l/iOCGlrFCUpghO7S3OZdlsha9E+yOkptaBbTBouKUcxzv6zM9Zi/9g5u2b3r+ZqqhyBlwatelsrnWL0x0jm7l1VlVYflUfVm559lDCnvBrk3VFM8ArR1tf3f7cy8+nx5gQ7YTfk/y+7Yhf0Aqt+K07XZfpIZX/UNjGFrLrLpXlK/UeZRkO548Rs4CvyKGIFi1hGGqvErYSeiwoJ0tpIKid7wZ9TiVPPC/oS5KaTqS/JmfXp+MJSWOxefMYYJqmis1jAN6yRN0FOf+EYf7SM1tAOtMN2qNYRaj7hSdN+PgBhkJ/M</latexit><latexit sha1_base64="bjxjqj8nvVAIAbZfa8fXZ/sgSRY=">AAADNHichVFNaxNBGH6y9aPGj0a9CF4GQyVRSCcVUQSlWgQvxTZt2kC3LLObSTJkv9zdBOOSP+Af8OBJwYN49w948S6CBW/iQQQvFUTw4LuTlaqlOsPMvPPM+zzzzLx26Ko44Xy7YEwdOHjo8PSR4tFjx0/MlE6eWo+DQeTIphO4QdSyRSxd5ctmohJXtsJICs925YbdX8zON4YyilXgryWjUG55ouurjnJEQpBV8phpy67yU3lvoKELY1ZkZthTrGIKN+yJ6vVOxbQDtx2PPFpSM+nJRIytPrvIJhkk4bWtftVk5i8Su8F4lXSk3/5d2SqVeY3rxvYG9TwoI2/LQekdTLQRwMEAHiR8JBS7EIipb6IOjpCwLaSERRQpfS4xRpG4A8qSlCEI7dPcpd1mjvq0zzRjzXboFpdGREyGWf6WP+c7/DV/wT/xH/tqpVoj8zKi1Z5wZWjNPDyz+u2/LI/WBL1d1j89J+jgqvaqyHuokewVzoQ/fPBoZ/VaYzY9z5/yz+T/Cd/mr+gF/vCr82xFNh5r9V1uijk09V0RzXP6P7o0RtqdIGQNt9GChSXcxCIauEvYUu7Z0pFD3AA17eR+9udU4vrfBd0brM/X6pdql1fmywu38mJP4yzOoUIVvYIF3MEyeXPwBt8LRmHKeGm8Nz4YHyepRiHnnMYfzfjyE1o1xEU=</latexit>を考える．一変数関数の最適化問題を考える(0) = f (✓ k )<latexit sha1_base64="PwFPMy2QBC+91JN5gX2KyLJnL24=">AAAC73ichVHPaxNBFP66/qqp2qgXwUswqaSX9KUiiiBUi+Cl2KZNG2jKsjudJEP2l7uTYFz6D3jw2oMgKngoXr3qpRf/AcH+CeKxghcPvp0sFC3qG2bem2/e9+abeW7kqUQTHUxYJ06eOn1m8mxh6tz5C9PFi5fWk3AQC9kUoRfGLddJpKcC2dRKe7IVxdLxXU9uuP3F7HxjKONEhcGaHkVyy3e6geoo4WiG7GK1Uim0o56q0mzpbqlTbbuht52MfHZpW/ekdnbs/myhUrGLZaqRsdLxoJ4HZeS2HBa/oI1thBAYwIdEAM2xBwcJj03UQYgY20LKWMyRMucSOygwd8BZkjMcRvu8dnm3maMB77OaiWELvsXjGTOzhBn6THt0SJ/oHX2ln3+tlZoamZYRe3fMlZE9/ezK6o//snz2Gr0j1j81a3Rw22hVrD0ySPYKMeYPn+4ert5pzKTX6Q19Y/2v6YD2+QXB8Lt4uyIbL0z1I26KOTTNXTGvc+Y/ujxHRp3DyBoeoAUbS7iHRTTwiLGlXLNtIsHcEDWj5En259zi+p8NPR6sz9fqN2o3V+bLC/fzZk/iKq6hyh29hQU8xDJrE3iO9/iAj9Zja9d6ab0ap1oTOecyfjNr7xcfdKrM</latexit>

よりも目的関数の値が小さい(0) > (↵)
<latexit sha1_base64="oA0Az0xkxc+1RTlbguVuJ/U7ryA=">AAAC43ichVFBTxNBGH2sglARil5M5NDYYuBSvmKIxoMBCYkXAhQKTSjZzI5DO2G7u9ndNtaGC0e9GRMPxoMmJhpO+he8+AdM5CcQjpBw8eC3002IEnUmM/Pmzfe+eTOfE7g6iokO+6xLl/sHrgwOZa4OXxsZzY5d34j8VihVRfquH1YdESlXe6oS69hV1SBUoum4atPZXUjON9sqjLTvrcedQG03Rd3TO1qKmCk7O14oZGpBQ0/SVO5hzqCacIOGmMoUCnY2T0UyLXcRlFKQR9pW/OwP1PAEPiRaaELBQ8zYhUDEfQslEALmttFlLmSkzbnCHjKsbXGU4gjB7C7Pdd5tpazH+yRnZNSSb3F5hKzMYYK+0yc6oW90QEf086+5uiZH4qXDq9PTqsAefX5z7ey/qiavMRrnqn96jrGD+8arZu+BYZJXyJ6+/ez1ydqD8kT3Dr2nY/b/jg7pK7/Aa5/KD6uq/MZkP9d2MY2KuSvkedr8R51Hx7gTzKxjEVXYWMI8FlDGMnNLqWfbIMlaH0Xj5Gny51zi0p8FvQg2Zoqlu8XZ1Zn83KO02IO4hduY5IrewxweY4W9SezjIz7ji6WsF9ZL61Uv1OpLNTfwW7Pe/gJ5J6Un</latexit>となる↵ を上手く選ぶ<latexit sha1_base64="lR1tPdba67cuLJznjCrGeU9dE7s=">AAAC0XichVFNSxtBGH5cP2rjR6K9CL2IieIpTiyiePKDQi+iJiYGjITZcZIsbnaX3U0whoJ4bO/twYOt4EH8GV76Bwr6E8RjBC8efHeyICrqDDPzzDPv884z8+qOaXg+Y1cdWmdXd8+H3o+Rvv6BwWhsaDjn2TVXyKywTdvN69yTpmHJrG/4psw7ruRV3ZSb+u5ycL5Zl65n2NaG33DkdpWXLaNkCO4TVUgkIgVuOhUeSSSKsThLMtVGX4JUCOII25odu0QBO7AhUEMVEhZ8wiY4POpbSIHBIW4bTeJcQoY6l/iOCGlrFCUpghO7S3OZdlsha9E+yOkptaBbTBouKUcxzv6zM9Zi/9g5u2b3r+ZqqhyBlwatelsrnWL0x0jm7l1VlVYflUfVm559lDCnvBrk3VFM8ArR1tf3f7cy8+nx5gQ7YTfk/y+7Yhf0Aqt+K07XZfpIZX/UNjGFrLrLpXlK/UeZRkO548Rs4CvyKGIFi1hGGqvErYSeiwoJ0tpIKid7wZ9TiVPPC/oS5KaTqS/JmfXp+MJSWOxefMYYJqmis1jAN6yRN0FOf+EYf7SM1tAOtMN2qNYRaj7hSdN+PgBhkJ/M</latexit>

二分法のアルゴリズムとなる ↵1 を定める．

1． ↵0 = 0 とし，<latexit sha1_base64="bApT8KexhyVAQGttxRqWUb/cVwQ=">AAAC1nichVFNSxtBGH6yflTjR6JeBC/BxOIpzipF8aSVgpdgEhMNqCyz6xgXN7vL7iYYQ3srhZ5685CTogfxZ3jpHyjoTxCPFnrx4LuTBamizjAzzzzzPu88M6/uWqYfMHYTU7q6e3o/9PXHBwaHhhPJkdEN36l7higbjuV4FZ37wjJtUQ7MwBIV1xO8pltiUz9YCc83G8LzTccuBU1X7NR41Tb3TIMHRGnJRCYT3+aWu881NZ7JpLRkmmWZbKmXQI1AGlHLO8lrbGMXDgzUUYOAjYCwBQ6f+hZUMLjE7aBFnEfIlOcCXxEnbZ2iBEVwYg9ortJuK2Jt2oc5fak26BaLhkfKFKbYH3bB7tlvdslu2cOruVoyR+ilSave0QpXS/wcX//3rqpGa4D9J9WbngPsYUF6Ncm7K5nwFUZH3zg6vl9fLE61PrJTdkf+T9gNu6IX2I2/xnlBFNsy+5O2hRmU5V0ezTPyP6o0mtIdJ6aEL6hAQw7LWEERa8TlIs+aRAZpHWSlk8Pwz6nE6vOCvgQbs1l1LvupMJte+hwVuw8TmMQ0VXQeS1hFnryFVW/jDOdKRfmmfFd+dEKVWKQZw39N+fUIFE+gyw==</latexit>とする．

<latexit sha1_base64="bu9zDGavkIX9I6hhpgU8sSP7crA=">AAAC2nichVHPSxtBGH2u1h+p1agXoZeliaWn+MUiiiBYpdCLVKPRgJFldjrGxc3usrsJpsGLN/Va6KGHUqFQ6Z/Ri/+AUP8E8WihFw/9drIgrbSdYWbevPneN2/mswPXiWKiyy6ju+dBb1//QObh4KOh4ezI6EbkN0KpytJ3/bBii0i5jqfKsRO7qhKEStRtV23ae0vJ+WZThZHje+txK1DbdVHznB1HipgpKzuWz2eqwg12hUXmvEmZfN60sjkqkG7mfVBMQQ5pW/Gz31HFG/iQaKAOBQ8xYxcCEfctFEEImNtGm7mQkaPPFQ6QYW2DoxRHCGb3eK7xbitlPd4nOSOtlnyLyyNkpYkJuqAzuqFz+kpXdPvXXG2dI/HS4tXuaFVgDR+Pr/38r6rOa4zdO9U/PcfYwaz26rD3QDPJK2RH33z7/mZtrjTRfkqndM3+P9ElfeMXeM0f8vOqKn3Q2e+0bUyirO8KeZ7U/1Hj0dLuBDPreIkKLCzjBZZQwmvmllPPlkaStT4K2sl+8udc4uKfBb0PNqYKxeeF6dWp3MJiWux+PMYTPOOKzmABr7DC3iR7+YgvODOqxqFxZJx0Qo2uVDOG35rx7hdW+6Gf</latexit>↵k 2 + ↵k 1
2． ↵k =
2とする

<latexit sha1_base64="xmwD38Gn3WOgH2eQOJc/U9Lesko=">AAAC/nichVHPaxNBFP662h9GbaNeBC9Lk4pQTGcjUhGE2iL0UmzTpg00ZZmdTtIlm91ldxOMy4Jn/4EexIOCBZHi1XsvHvQo2D9BPLbgxYNvJwuhLbVvmJn3vnnfm2/mWb5jhxFjh0PapcvDI6NjV3JXr10fn8jfuLkeep1AyKrwHC+oWTyUju3KamRHjqz5geRty5EbVmshPd/oyiC0PXct6vlyq82brt2wBY8IMvOzxWKuzh1/h5st/YlebwRcxBkQt+6XE31aH4RGksTlJFcs6ma+wEpMmX7WMTKngMyWvfxP1LENDwIdtCHhIiLfAUdIYxMGGHzCthATFpBnq3OJBDnidihLUgYntEVrk6LNDHUpTmuGii3oFodmQEwdU+wH+8iO2Ff2if1if8+tFasaqZYe7VafK31z4vXt1T8Xstq0R9gZsP6rOUIDj5RWm7T7CklfIfr87svdo9XHlan4LnvPfpP+d+yQHdAL3O6x+LAiK29U9QE3xgyq6q6A1hn1H02aPaWOE7KGZ6jBxBKeYgEVPCdsKdNsKk8Q10NJKXmR/jm12Djd0LPOerlkPCg9XCkX5uazZo/hDiZxjzo6izksYpm0CbzFAb7hu/ZK29P2tc/9VG0o49zCCdO+/APROLCc</latexit>

3．停止条件が満たされるならば，↵k を数値解として出力して停止<latexit sha1_base64="D3CJKZCkjVJqC3Rgs/hOqGX8Owo=">AAAC1nichVHPSxtBGH2utdVUm6gXwUtoEukpmSjF0pM/ELxINTEaSMIyu07iks3usrsJTYPeRPDkzYMnRQ/in9GL/4DQ/AniUcFLD/12siA22M4wM2/efO+bN/Npjml4PmPdAWXwzdDbd8MjkfejYx+isfGJLc9uuroo6LZpu0WNe8I0LFHwDd8URccVvKGZYlurLwfn2y3heoZtbfptR1QavGYZVUPnPlFqLJpMRsrcdHa5Wo8kk3E1lmBpJlu8H2RDkEDY1u3YL5SxAxs6mmhAwIJP2ASHR72ELBgc4iroEOcSMuS5wB4ipG1SlKAITmyd5hrtSiFr0T7I6Um1TreYNFxSxpFit+yKPbAbds3u2O9Xc3VkjsBLm1atpxWOGj2ayj/9V9Wg1cfus+qfnn1U8UV6Nci7I5ngFXpP3/px8pD/mkt1Ztg5uyf/Z6zLftILrNajfrkhcqcy+7O2gwwK8i6X5oz8jxqNtnTHidnECopQsYZFLCOHb8SthZ5ViXTS2khLJ9+DP6cSZ/8uaD/Ymk1n59KfN2YTC0thsYcxjY/4RBWdxwJWsU7egqqf4gKXSlHZVw6Uw16oMhBqJvGiKcd/AJ9roQU=</latexit>

04．
<latexit sha1_base64="JqUMe7yhF/RmMEaGA3uYvBlR1WU=">AAADB3ichVFNaxNBGH66ftX40dReBC+DSbUemk5SxCIK0SJ4KbZp0waassyu02TIZneZ3UTj0h+gf8CDJwUF8eBNvIoXD14r9ieIxwpePPjuZEvRop1hZp553nneeWZeJ/RUFHO+M2IdOXrs+InRk7lTp8+cHcuPn1uNgp52Zd0NvEA3HBFJT/myHqvYk41QS9F1PLnmdObT+Fpf6kgF/ko8COVGV7R8talcERNl56vFYo41w7a6zKaawgvbwu5cYTcYZ02tWu1YaB08YHsRdnMPJp3pylauWGR2vsBL3DR2EJQzUEDWFoP8VzRxHwFc9NCFhI+YsAeBiPo6yuAIidtAQpwmpExcYgs50vbolKQTgtgOzS3arWesT/s0Z2TULt3i0dCkZJjk2/w13+Wf+Bv+jf/6Z67E5Ei9DGh1hloZ2mNPzi//PFTVpTVGe1/1X88xNjFnvCryHhomfYU71PcfPd1dvl6bTC7xF/w7+X/Od/hHeoHf/+G+XJK1Zyb7vjbBDOrmLk3zjPmPFo2BcSeIWcEdNGBjAbcwjxruEbeQebYNckkboGScPEz/nEpc/rugB8FqpVSeLV1dqhSqt7Nij+ICLmKKKnoNVdzFInlz8QqfsY0v1mPrrfXOej88ao1kmgn80awPvwGIM7MB</latexit><latexit sha1_base64="siLE14p53iI9fabHFR34P8b1QEw=">AAADB3ichVFNaxNBGH669aPGj0a9CF4Gk2o9NJ2tSEVQokXwUmzTpg00ZZldp8mQze4yu4nGpT9A/4AHTwoK4sGbeBUvHrxW7E8QjxW8ePDdyZaiRTvDzDzzvPO888y8buSrOOF8e8QaPXT4yNGxY4XjJ06eGi+ePrMShz3tyboX+qFuuCKWvgpkPVGJLxuRlqLr+nLV7cxl8dW+1LEKg+VkEMn1rmgFakN5IiHKKVbL5QJrRm11iU02hR+1hdO5zG4yzppatdqJ0Dp8wHYj7MYuTDtT9mahXGZOscQr3DS2H9g5KCFvC2HxK5q4jxAeeuhCIkBC2IdATH0NNjgi4taREqcJKROX2ESBtD06JemEILZDc4t2azkb0D7LGRu1R7f4NDQpGSb4Fn/Nd/gn/oZ/47/+mSs1OTIvA1rdoVZGzviTc0s/D1R1aU3Q3lP913OCDVwzXhV5jwyTvcIb6vuPnu4sXa9NpBf5C/6d/D/n2/wjvSDo//BeLsraM5N9T5tiGnVzl6Z52vxHi8bAuBPELOMOGnAwj1uYQw33iJvPPTsGeaQNUTFOHmZ/TiW2/y7ofrAyU7GvVK4uzpSqt/Nij+E8LmCSKjqLKu5igbx5eIXP2MIX67H11npnvR8etUZyzVn80awPvwGK5LMC</latexit>(↵k ) < 0 ! ↵k = ↵k 20(↵k ) > 0 ! ↵k = ↵k 1

0<latexit sha1_base64="X0sGEJ6LVVWJHU4VdAjBdeagUbY=">AAAC+HichVHPaxNBFP66/mrjj6Z6EbwMJtV6SWcrRRGE1iJ4KbZp0waassyO02TIZneZ3UTT0H+g/0APXlQQKd68exLUf0Cwf4J4rODFg28nK0WL+oaZefO99735Zp4fBzpJOT8YcU6cPHX6zOhY4ey58xfGixMX15Koa6SqySiITN0XiQp0qGqpTgNVj40SHT9Q6357IYuv95RJdBSupv1YbXZEM9RbWoqUIK/olssF1ohb+jqbaoggbgmvfYPdZZw1jG62UmFM9Jj9ihTKZeYVS7zCrbHjjps7JeS2FBU/o4FHiCDRRQcKIVLyAwgkNDbggiMmbBMDwgx52sYVdlAgbpeyFGUIQtu0Num0kaMhnbOaiWVLuiWgaYjJMMk/8X1+yD/y1/wL//HXWgNbI9PSp90fclXsje9eXvn+X1aH9hStI9Y/NafYwm2rVZP22CLZK+SQ39veO1y5U50cXOMv+FfS/5wf8Hf0grD3Tb5cVtWntvoRd4Bp1OxdhtZp+x9Nmn2rThCyivuow8Mi5rGAKh4Stphr9qwniRuhYpU8yf6cWuz+2dDjztpMxb1ZmV2eKc3dy5s9iiu4iinq6C3M4QGWSJvEHt7iPT44284z55WzP0x1RnLOJfxmzpufUqmtbA==</latexit>5．

(↵k ) = 0 ! ↵kを数値解として出力して停止として，ステップ2に戻るサンプルプログラムのダウンロード

パラメータ推定のサンプルプログラムは以下からダウンロードしてください

http://mmc01.es.hokudai.ac.jp/els/link
数理モデルに対するパラメータ推定の解説はこちらをご覧下さい．
サンプルプログラムはこちらからダウンロードしてください．

サンプルプログラム#ロジスティック方程式 最急降下法python
import numpy as np
import matplotlib.pyplot as plt
import sys
import random#微分方程式の右辺
def f(x,par):return par[1]*(par[2]-x)*x#Runge-Kutta法
def runge(x,par):k1 = dt*f(x,par)k2 = dt*f(x+0.5*k1,par)k3 = dt*f(x+0.5*k2,par)k4 = dt*f(x+k3,par)x += (k1+2*k2+2*k3+k4)/6
return x#コスト関数def costfunc(par,y,st):xc=[]
x=par[0]for i in range(N):x=runge(x,par)#時刻情報stと同じiのループの時にシミュレーションデータとして保存
if i in st:
xc.append(x)s=0.0
for i in range(Nz):s += (xc[i]-y[i])**2
return(s+r*sum([x**2 for x in par]))サンプルプログラム#コスト関数の勾配を計算
def grad(par,y,st):par̲grad=[]for m in range(len(par)):par̲plus=[par[m]+0.01*par[m] if i==m else par[i] for i in range(len(par))]par̲minus=[par[m]-0.01*par[m] if i==m else par[i] for i in range(len(par))]par̲grad.append((costfunc(par̲plus,y,st)-costfunc(par̲minus,y,st))/(2*par[m]/100.0))
return par̲grad

サンプルプログラム#二分法(直線探索)
def bisection(par,par̲grad,y,st):##初期設定s=0.0#二分法左端初期値
l=0.01#二分法右端初期値
count=0g̲plus=[par[i]-(l+0.01*l)*par̲grad[i] for i in range(len(par))]g̲minus=[par[i]-(l-0.01*l)*par̲grad[i] for i in range(len(par))]phiplus=costfunc(g̲plus,y,st)phiminus=costfunc(g̲minus,y,st)#計算可能かつ現在のコスト関数より小さいlになるまでlを小さくする
while np.isinf(phiplus) or costfunc(par,y,st)<phiminus :
l=l/2.0g̲plus=[par[i]-(l+0.01*l)*par̲grad[i] for i in range(len(par))]g̲minus=[par[i]-(l-0.01*l)*par̲grad[i] for i in range(len(par))]phiplus=costfunc(g̲plus,y,st)phiminus=costfunc(g̲minus,y,st)count += 1
if count>100:
print("ERROR̲bisection")sys.exit()#右側の微分が負になるまでlを右にずらしていく
while (phiplus-phiminus)<0.0:
old̲l=ll=1.1*lg̲plus=[par[i]-(l+0.01*l)*par̲grad[i] for i in range(len(par))]g̲minus=[par[i]-(l-0.01*l)*par̲grad[i] for i in range(len(par))]phiplus=costfunc(g̲plus,y,st)phiminus=costfunc(g̲minus,y,st)サンプルプログラム#計算が発散したら一つ前のlを出力
if np.isinf(phiplus) or np.isinf(phiminus):
print("ERROR̲l=%.15e" %l)return old̲lwhile True:
p=(s+l)/2.0#二分法の右端と左端の真ん中の点をpとする
g̲plus=[par[i]-(p+0.01*p)*par̲grad[i] for i in range(len(par))]g̲minus=[par[i]-(p-0.01*p)*par̲grad[i] for i in range(len(par))]phiplus=costfunc(g̲plus,y,st)phiminus=costfunc(g̲minus,y,st)if (phiplus-phiminus)<0.0:#pでの傾きが負だったらpを左端にする
s=pelif (phiplus-phiminus)>0.0:#pでの傾きが正だったらpを右端にする
l=pelif (phiplus-phiminus)==0.0:
l=s=p
else:
print("ERROR!!")
sys.exit()
if(abs(l-s)<1e-8):#区間が十分小さくなったら終了
p=(s+l)/2.0
return p
print("ERROR! bisection¥n")sys.exit()

サンプルメイン関数################################################メイン文開始#################################################変更しないパラメータ
dt=0.01 #刻み幅
N=500 #ルンゲクッタ試行回数
r=0.0000001 #正則化パラメータ#自分で実験データを作る場合 st,yを使用######################################
st=[25,100,250,470] #データ点の時刻情報
Nz=len(st)#データ点の個数取得
par̲ans=[0.3,1.0,2.0] #正解パラメータ，左からx0,a,b#正解の解軌道t=0.0
tpoints=np.arange(0.0,(N+1)*dt,dt)
xpoints=[]x=par̲ans[0]y=[]xa=[] #ノイズ用
for i in range(N+1):
xpoints.append(x)x=runge(x,par̲ans)
xa.append(random.uniform(0.95*x,1.05*x)) #ノイズ用
if i in st:y.append(x) ← ここでy.append(xa[i] ) とすると10％の誤差が混入できる

#実験データyの表示
map̲y=map(str,y)y̲show=', '.join(map̲y)print("y="+ y̲show)###########################################最急降下法####パラメータ初期値
par=[1.0,3.0,5.0]#カウント数n̲count=0
while True:
if n̲count>10000:
print("looperror")sys.exit()
z=costfunc(par,y,st) #現在のコスト関数の値
print("n=%d costfunc=%.15e" %(n̲count,z))par̲grad= grad(par,y,st) #グラディエント計算
alpha = bisection(par,par̲grad,y,st) #直線探索ステップ幅計算
par̲old=par #一つ前のパラメータを保存
par=[par[i]-alpha*par̲grad[i] for i in range(len(par))] #パラメータ更新
s=0.0
for i in range(len(par)):s += (par[i]-par̲old[i])**2
if s<1e-8: #パラメータの更新が終了判定を満たしたら終了#計算にかかった更新の回数 コスト関数 導出されたパラメータ
print("n=%d" %n̲count)
print("costfunc=%f" %costfunc(par,y,st))print(par)breakn̲count += 1####主な計算終了

#導出されたパラメータの解軌道
xpoints̲a=[]x=par[0]for i in range(N+1):
xpoints̲a.append(x)x=runge(x,par)#グラフ描写plt.plot(tpoints,xpoints)
plt.plot(tpoints,xpoints̲a)for i in st:
plt.plot(i*dt,xpoints[i],marker='o',color='r',ms=10)# plt.plot(i*dt,xa[i],marker=ʻoʼ,color=ʻrʼ,ms=10) #ノイズ用
plt.show()最適化問題の結果
取得定量データの個数による誤差変化
例：ロジスティック方程式

N=3

N=5

N = 10N = 15

パラメータは取得データ数を増やす事で，パラメータ推定の精度は上がることが分かった．
また，データにノイズを加えた場合，取得データの個数が少ないときは解軌道
はノイズを与えたデータ上をなぞるように動くが，取得データを増やすことでノイズを与えたデータの影響が小さくなり，与えたモデルの解軌道に近づくことが分かった.

データ点の位置による影響
取得定量データの時間位置による誤差変化

A

B

C

データを同じ個数取ってきても,データの取得位置によって解軌道に部分的な差が生じる.

D

データを取得してくる時間位置毎に,取得データがどのパラメータの情報を多く有しているかが異なる.

サンプルプログラムを動かそう！
ロジスティック方程式に対するパラメータ推定

プログラム名: sample01.py
ターミナルを立ち上げて，コマンドラインから% python3 sample01.pyとすれば計算でき，結果がグラフとして出力されます．

パッケージ不足のエラーが出たら以下のコマンドからパッケージをインストールする% pip3 install パッケージ名

それでも動かなかったら…，TAに聞きましょう！

演習問題
2種競争系モデル

初期値

数理モデルの設定から
u0, v0, b, c は非負の実数データから決定したいパラメータ

sample01.pyを参考にしてプログラムを書いてみましょう！

答え% python3 sample02.py

1．

のとき，単安定，例えば，b=2.0, c=0.52．

のとき双安定．b=0.5, c= 2.0

ベイズ推定を用いたパラメータ推定
データ点の取り方によって，モデル方程式から生成したデータにも関わらず最適化問題では収束してもパラメータを確定できない場合がある．
（ある意味データの性質が悪いとも言える）→パラメータ同定不能問題→谷口さん（神戸大学）の研究：グレブナー基底を使った解析からパラメータ同定不能条件を示している
パラメータ同定が不可能な場合を込めてベイズ推定を使ってパラメータを確率分布として求めてみる．

ベイズ推定：機械学習やデータ解析等でよく用いられる方法

ベイズの定理
を定数， を確率変数とする．ベイズの定理より以下の関係式が成り立つ．

ベイズ推定とは事前分布(結果)と尤度関数からベイズの定理を用いて事後分布(原因)を推定することである．事後確率を確率分布として求めることでパラメータ推定を行う．

ベイズ推定の適用
データの個数を，数理モデルのパラメータを ,実験データをパラメータ,
下でのシミュレーションデータをとする．

に対する との相対誤差は正規分布：
以下の式が成り立つと仮定する
ここで，

は分散

，平均0の正規分布に従う確率変数.

このとき，尤度関数は以下のように定義できる：

実験データとシミュレーションデータの距離が近いほど良いという指標(尤度が大きい)に従うとし，

ベイズの定理再出

ベイズの定理

事後分布

尤度関数

事前分布

推定
事前分布にはパラメータに関してのアプリオリな情報を取り込む．
事前に得ている情報が少ない時は，広い区間での一様分布や分散がとても大きいガウス分布を用いる．(無情報事前分布)

事後分布「与えられた実験データに対して数理モデルが尤もらしい解を生成するパラメータの確率」

事後分布の推定方法
マルコフ連鎖モンテカルロ法(MCMC)

MCMC

多変量の確率分布からの乱数発生アルゴリズム

メトロポリスヘイスティング法（M-H法）

MCMCの代表的な手法の一つ

M-H法を用いて事後分布からのサンプリングを行う

M-H法アルゴリズム（ベイズ推定）

求めたい分布1. 初期値からのサンプリングを決める．
に対して次を繰り返す．

2.

を提案分布
から発生させる.(提案分布はの一様分布から発生させを区間を中心とした一様分布)の場合そのほかの場合

遷移確率とする．
遷移確率の計算(

提案分布の対称性)(

ベイズの定理)M-H法アルゴリズム（ベイズ推定）改訂
求めたい分布1. 初期値のサンプリングを行う
を決める．

2.

に対して次を繰り返す
を提案分布（提案分布はを区間から発生させる．
を中心とした一様分布)の一様分布から発生させる．

遷移確率

とする．
遷移確率の計算(

提案分布の対称性)(

ベイズの定理)M-H法の利点
最急降下法→

次の候補点をその点の勾配によって決定
局所最適解に陥りやすい

M-H法

→

次の候補点をランダムウォークによって決定
局所最適解に陥りにくい

M-H法によるサンプリングのイメージ

ベイズ推定を用いた数理モデルのパラメータ推定の適用例

パラメータをとする．
下図の青い点を実験データとし，これらの点を用いて実験1,2の場合におけるを推定する．

ロジスティック方程式

赤線：解軌道 青点：実験データ
2

2

1.8

1.8

1.6

1.6

1.4

1.4

1.2

1.2

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

実験1：初期状態から平衡状態までほぼ均等な位置の実験データ6点

0.2

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

実験２：平衡解付近のみの実験データ4点

ベイズ推定を用いた数理モデルのパラメータ推定の適用例
1.4

3

1.2

2.5

1

実験1

0.6
0.4

3

1.5

2.5

1

1.5

2
1

0.5

0.2
0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

0

0.5
1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5

0.11

3.5

0.1

3

0.09

2.5

0.08

実験2

4
3.5

2

0.8

0

4.5

0.07

2

0.06

1.5

0.05

1

0.04

0.5

0.03
0.02
-1

0

1

2

3

4

5

6

7

8

9

10 11

0
-0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5

0
-0.3
0.12
0.11
0.1
0.09
0.08
0.07
0.06
0.05
0.04
0.03
0.02
-1

0

0

1

0.3

2

3

0.6

4

5

0.9

6

7

8

(各パラメータのサンプリングをヒストグラムにしたもの)

パラメータの真値
サンプリングは全て
メタパラメータ は個，事前分布はとした．

の一様分布を仮定し，

1.2

1.5

9

10 11

パラメータ間の相関図

実験1におけるサンプリングされたパラメータ間の相関係数のヒートマップと散布図

実験2におけるサンプリングされたパラメータ間の相関係数のヒートマップと散布図

相関係数
データ列

ただし

と

の標本相関係数サンプルプログラム
MHMCを用いたベイス推定の計算（ロジスティック方程式）% python3 sample02.py
パラメータの分布を可視化% pythos3 makehistogram.py
相関係数とデータの分布図を計算% cat xxxxx.dat | python3 heatmap.py

演習問題
2種競争系モデル

初期値

数理モデルの設定から
u0, v0, b, c は非負の実数データから決定したいパラメータ

sample02.pyを参考にしてプログラムを書いてみましょう！

答え% python3 mcmc̲Competition.py1．

のとき，単安定，例えば，b=2.0, c=0.52．

のとき双安定．b=0.5, c= 2.0

参考文献

金森敬文，鈴木大慈，竹内一郎，佐藤一誠著，
機械学習のための連続最適化，講談社

C.M.ビショップ著，パターン認識と機械学習下，丸善出版

数理モデルの数値計算法について

数値計算に必要な知識は以下から！

http://mmc01.es.hokudai.ac.jp/els/link「私にとっての反応拡散系数値シミュレーション入門，基礎編」のPDFファイルのダウンロード．