今回は乱入(?)ネタ。9月に井上さんの書いた IGOR Pro の微分方程式ネタの記事についたコメントにあった、「微分方程式中のパラメーターをフィッティングで求められないか」という話題に対する Mathematica 的な解法の紹介。
もっと早く乱入すべきだったんだけど、最近しばらく心を落ち着けて検討することができなかったので、ようやく着手した次第。ちなみに以下は、最近日本語版がリリースされた ver.8 で行なった結果だが、ver.7 でも同じ方法で実行可能だ。

モデルの系は、井上さん記事で紹介されたのと同じものを使おう。 Mathematica で連立微分方程式モデルは、割と見たままの表記で入力することができる。こんな感じ。

In[1]:= model = {a'[t] == b'[t] == -k1 a[t] b[t] + k2 c[t],
          c'[t] == k1 a[t] b[t] - k2 c[t] - k3 c[t],
          d'[t] == k3 c[t],
          a[0] == a0, b[0] == b0, c[0] == c0, d[0] == d0};

組み込み以外のシンボルは小文字から始めるという Mathematica の通常流儀にならって A 〜 D ではなく a 〜 d を使っていること、a 〜 d を t の関数の形で書いていること、t の微分は ' で書けること、a[0] 〜 d[0] の初期条件値も a0 〜 d0 というパラメーターでモデルに含めていること、くらいがポイントだろうか。

フィッティングの話の前に、まず IGOR Pro と同じように、パラメーターが決まっているものを数値的に解く方法から。
初期条件値として a[0] = b[0] = 1, c[0] = d[0] = 0 を、パラメーターは k1 = 0.02, k2 = 0.01, k3 = 0.04 を使い、時間範囲 0 〜 500 の数値解を得るには、微分方程式を数値的に解いてくれる NDSolve を使って次のようにする。

In[2]:= cons = {a0 -> 1, b0 -> 1, c0 -> 0, d0 -> 0};
        sol = NDSolve[
          model /. cons /. {k1 -> 0.02, k2 -> 0.01, k3 -> 0.04},
            {a, b, c, d}, {t, 0, 500}]

Out[3]= {{a -> InterpolatingFunction[{{0., 500.}}, <>],
          b -> InterpolatingFunction[{{0., 500.}}, <>],
          c -> InterpolatingFunction[{{0., 500.}}, <>],
          d -> InterpolatingFunction[{{0., 500.}}, <>]}}

結果はこのような補間関数 InterpolatingFunction の組の形で得られる。今回の系では1組しか解がないが、系によっては解が複数組得られることがあるため、結果がリストのリストの形になっていることに注意。

得られた結果をグラフにしてみよう。今回の系では a[t] と b[t] は同一の曲線になるので、a[t], c[t], d[t] を描いてみる。IGOR Pro 同様、中間体 c[t] が最初増えてから減っていく、化学の教科書に出てきそうなグラフを描くことができる。

In[4]:= Plot[Evaluate[{a[t], c[t], d[t]} /. First[sol]], {t, 0, 500}]

Out[4]= Mathematica出力画像 Out[4]

さて、今回の本題である、「微分方程式の中のパラメーターをフィッティングで求める話」に移ろう。
「フィッティングで」ということは、なんらかの観測データがあるということだ。ここでは、最終生成物の時間変化 d[t] の観測データがあるという前提で、次のようなデータを用意してみた。

In[5]:= data = {{0., 0.001}, {20., 0.087}, {40., 0.236}, {60., 0.357},
          {80., 0.456}, {100., 0.527}, {120., 0.584}, {140., 0.629},
          {160., 0.672}, {180., 0.697}, {200., 0.722}, {220., 0.740},
          {240., 0.760}, {260., 0.782}, {280., 0.794}, {300., 0.810},
          {320., 0.815}, {340., 0.823}, {360., 0.828}, {380., 0.837},
          {400., 0.859}, {420., 0.853}, {440., 0.866}, {460., 0.869},
          {480., 0.880}, {500., 0.872}};
        dplot = ListPlot[data]

Out[6]= Mathematica出力画像 Out[6]

フィッティングするためには、フィッティング対象になる関数が必要だ。これを、フィッティングパラメーターと t の関数 myF として定義する。上で書いたのと同じように、時間範囲 0 〜 500 で NDSolve した結果から d[t] の部分を取り出す内容だ。

In[7]:= Clear[myF];
        myF[inT_, inK1_?NumericQ, inK2_, inK3_] :=
          (d /. First[NDSolve[
            model /. cons /. {k1 -> inK1, k2 -> inK2, k3 -> inK3},
              {a, b, c, d}, {t, 0, 500}]])[inT]

実際のフィッティングは FindFit 関数で。 NonlinearModelFit 関数を使う手もあるが、 NonlinearModelFit はさまざまな統計テーブルを準備するためか時間がかかるので、こういう用途では FindFit がおすすめ。

In[9]:= myfit = FindFit[data, myF[t, k1, k2, k3],
          {{k1, 0.02}, {k2, 0.01}, {k3, 0.03}}, t]

Out[9]= {k1 -> 0.0182474, k2 -> 0.00758555, k3 -> 0.0426789}

うまくいくと、こんな感じで、フィッティング結果としてのパラメーターの変換規則のリストが得られる。

ポイントは2つ。
ひとつは、myF の引数定義で、フィッティングパラメーターのひとつに ?NumericQ をつけていること。これがないと、FindFit が記号的な導関数が得られないかを試そうとして(パラメーターに数値ではなくシンボルを与えて myF を評価して)、エラーが出まくった果てに失敗する。
もうひとつは、パラメーターの初期値。上の例では、探索する k1 〜 k3 の初期値として、0.02, 0.01, 0.03 を設定している。非線形最小二乗法では一般的に、パラメーター空間の探索初期値が極めて重要だ。現在検討している系にとって妥当な初期値、つまり、「データとばっちり一致するわけじゃないけど、だいたいこんなもんの値」という「それらしい値」を指定しないと、うまくいかないことは非常に多い。(探索初期値を省略すると、FindFit は 1.0 を初期値として採用する。それで上手くいくことももちろんあるのだが...)

得られた結果を使って、データと d[t] を重ね描きしてみる。少なくとも見た目の範囲ではうまくフィッティングできていることが分かる。

In[10]:= Show[Plot[Evaluate[myF[t, k1, k2, k3] /. myfit],
           {t, 0, 500}], dplot]

Out[10]= Mathematica出力画像 Out[10]

せっかく(?)なので、d[t] 以外の a[t] や c[t] とも重ね描きしてみよう。今回の myF の定義では、NDSolve の結果は保存されないので、a[t], c[t] を得るためには再度 NDSolve することになるが、こんな感じ。

In[11]:= sol2 = First[NDSolve[model /. cons /. myfit,
           {a, b, c, d}, {t, 0, 500}]];
         Show[Plot[Evaluate[{a[t], c[t], d[t]} /. sol2], {t, 0, 500}],
           dplot]

Out[12]= Mathematica出力画像 Out[12]

このグラフで赤色の線で描かれている中間体の濃度 c[t] が最大になる時間は、c'[t] がゼロ(c[t] の傾きがゼロ)になる時間。c'[t] は、最初に定義したモデルの第2式で与えられているので、これを使って次のように FindRoot すれば求められる。グラフから、探索開始時間として t = 20 を採用。

In[13]:= FindRoot[k1 a[t] b[t] - k2 c[t] - k3 c[t] /. myfit /. sol2,
           {t, 20}]

Out[13]= {t -> 27.811}

同様に、基質と生成物の濃度が逆転する時間(グラフの、青線と黄色線が交わる時間)も、a[t] == d[t] を解く形で次のようにすれば OK。

In[14]:= FindRoot[a[t] == d[t] /. myfit /. sol2, {t, 100}]

Out[14]= {t -> 79.171}

いかがだろうか。
今回は、NDSolve を使って解く連立微分方程式の系だったが、数値積分値を NIntegrate で求める関数のパラメーターをフィッティングするのも同じ方法で可能だ。また、独自に残差二乗和関数を定義して最小化する方法まで検討すれば、さらに多様な系のフィッティングを Mathematica で実行することが可能になる。

※この記事の内容は執筆者の個人的見解で、ヒューリンクスによる公式情報ではありません。[免責事項]

トラックバック

この記事へのトラックバックURL
http://blog.hulinks.co.jp/cgi/mt/mt-tb.cgi/452
内容に対しての関連性がみられないものは削除する場合があります

コメント一覧

先日「微分方程式中のパラメーターをフィッティングで求められないか」
という話題をコメントした張本人です。
たまに訪れておりましたが、思わぬお年玉でした。ありがとうございます。
やはり上記の場合はMathematica ですよね。
理解に欠けていた部分もあり、大変参考になりました。ありがとうございます。
上記の場合をegorで処理する事ができるなら非常に有難いとかねてから思っていました。
p.s.
井上様のegorからのアプローチも楽しみにしております。

守屋様の記事に対して返信する記事の場所を間違えてしまいました。
すいません。

重ねがさねコメント失礼します。
またお時間があれば
mathematicaによる色々なフィッティングも特集していただけると大変ありがたいです。
よろしくお願いいたします。

k2 さん、いつもありがとうございます。何かしらのご参考になったのでしたら、良かったです。
IGOR Pro 的な解法に関しては、近日中(たぶん)に、井上さんが新しい記事を掲載してくれるようですので、そちらをお待ち下さい。
Mathematica の追加フィッティングの話題については、Wolfram Research 社の下記ページなどに、さまざまな例が載っているので、参考になるのではないかと思います。

* FindFit の「アプリケーション」の項
http://reference.wolfram.com/mathematica/ref/FindFit.html#107405395
* NonlinearModelFit の「一般化と拡張」の項
http://reference.wolfram.com/mathematica/ref/NonlinearModelFit.html#1020263627

一歩踏み込んで「最適化問題」とすると、さらにたくさんのドキュメントが公開されています。例えば下記。
* 大域的非線形数値最適化
http://reference.wolfram.com/mathematica/tutorial/ConstrainedOptimizationGlobalNumerical.html

ともあれ、ありがとうございました。本年もよろしくお願いいたします。

この記事を参考にさせてもらっている某化学屋です。

こういった関数でのグローバルフィッテイングをする方法についても
よければ教えていただけないでしょうか?

NY さん、はじめまして。

グローバルフィッティングというと、複数のデータセットに対してパラメーター共有するような形でフィッティングするという話でしょうか。系に依存するのではないかと思います。
例えば、上記記事と似た形で初期濃度を変えるような系であれば、初期濃度の項を model にうまく持ち込むことができれば同じようにできるのではないかと想像します。model への組み込みが難しい場合は、FindFit を使わずに残差二乗和関数を定義して NMinimize を使うしかないのではないかと。

…なんかあんまり参考にならないコメントしかできなくてスミマセン。

便乗して質問させていただきたいのですが...
Findfitで得たパラメータ(本文中ではk1,k2,k3)を後に再利用したいと考えています。
具体的にはk1*k2などといった演算をして値を得たいと考えています。
しかし、単純にk1*k1などとしても新たに変数が定義されたと解釈されてしまうようです。
ご教示いただければ幸いです。

本記事を執筆した守谷は、なぜかヒューリンクスを巣立っていってしまったので、代わりに(わかる範囲で)お答えしますね。

FindFit が返すのは、本文にあるように変換規則のリストなので、フィッティングで得た値を使って k1*k2 を計算するには、/.(ReplaceAll)を使って k1*k2 を変換するのが普通だと思います。例えばこんな感じ。

k1*k2 /. myfit

シンボル k1, k2, k3 に FindFit の答えを直接代入してしまいたい場合は、こんな風にします。

{k1, k2, k3} = {k1, k2, k3} /. myfit

こうすれば、以後、k1, k2, k3 は、それぞれ、フィッテイングした値として使えます。
ただしこれをすると、同じ Mathematica セッション中に本文中の FindFit を再度実行するとエラーになるので(パラメーターがすでに値を持っている状態なので)、FindFit 再実行の前に、

Clear[k1, k2, k3]

をする必要がでてきます。

変換規則による置き換えについては、ヘルプをご参照ください。
http://reference.wolfram.com/mathematica/tutorial/ApplyingTransformationRules.html

井上様

場違いな質問にもかかわらずご回答下さいましてありがとうございました.

ご指摘のとおり記述したところ再利用することができました.

HDK さん

上手くパラメータが再利用できてよかったです。


コメントの投稿

Emailアドレスは表示されません。は必須項目です。
ヒューリンクス取り扱い製品の内容や購入に関するお問い合わせはヒューリンクスサイト連絡先へお願いいたします。投稿前にその他の注意事項もご覧ください。

HULINKS サイトの新着情報