多目的最適化について

質問

あるモデルが連立非線形微分方程式で表せるとします。
実際のデータとfittingすることによって、微分方程式の係数を、遺伝的アルゴリズムを用いて推定したいと思っています。

その時、目的関数(GAでの適応度)の指標として、今２つのものを考えています。（実際には、入力を加えるとパルス状の出力をするものなので、各点での誤差とパルスのタイミングの２つです）。

このような２つの量を目的関数に用いたい場合、どのような方法があるのでしょうか？？？
「重み付けして和を取るっていう方法」と「一方を制約条件として考える」以外で教えてください。
お願いします。

回答No.3

stomachman · ベストアンサー

●v(0)=0なんて嘘書いてしまいましたが、ちょっと混乱してました。忘れてください。
●面白そうだから計算をやってみようと思い、もう一度よく見たら、式の数値の単位が分からない。とりあえずmV, msと仮定して、グラフからv(t)のデータをえいやと起こし、m,h,nを計算してみたんですが、n(t)が全然似ない。何だか変ですよ。
an = {-0.01(V(t)+50)} /[exp{-(V(t)+50)+10}+1]
じゃなくて
an = {0.01(V(t)+50)} /[exp{-(V(t)+50)+10}+1]
では？
●ということにして、Δt =0.25msぐらい、というめちゃくちゃおおざっぱな計算をExcel(！)でやった。m0,h0,n0はHPに書いてあった数値を使いました。m(t)が多少振動しちゃったりしてますが、ま、それなりの答は出るようです。
具体的には線形最小二乗法
y=Ax+ε
ここにy,Aの縦ベクトルはそれぞれ
y(t) = integral{s=0～t} I(s)ds （～8E-3）
a[1](t) = v(t)
a[2](t) = integral{s=0～t} v(s)h(s)m(s)^3 ds
a[3](t) = integral{s=0～t} h(s)m(s)^3 ds
a[4](t) = integral{s=0～t} v(s)n(s)^4 ds
a[5](t) = integral{s=0～t} n(s)^4 ds
a[6](t) = integral{s=0～t} v(s)ds
a[7](t) = t
a[8](t) = 1
です。右辺はいい加減な台形則で積分。そして単純にεの二乗和を最小化するために
x = (A' A)~ y
ここに(A')は転置、~は逆行列です。
で、得られた数値は
C 0.000109906
Gnamax 0.002923773
-GnamaxEna -0.196207736
Gkmax 0.000448064
-GkmaxEk 0.037885549
Gl 0.000161473
-GlEl 0.010964163
constant 0.00986317
ちょっとおかしいけど、まあなんとなく計算できるらしいことは分かった。（最後のconstantってのは、方程式を積分しちゃったせいで出てくる余計な項の積もりですが、要らないのかも知れない。）あとは根性入れて細かくやれば旨く行くのかどうか。旨く行ったらヒトヲクッタ法てとこでしょうね。

お礼

遅くなりました。紹介していただいた方法も試してみたいと思います。
数値積分の刻み幅次第では比較的よい推定ができそうな気がしています。

式のことですが、紹介したホームページの式は
αm = {-0.1(V+35)}/[exp{-(V+35)+10}-1]
αn = {-0.01(V+50)}/[exp{-(V+50)+10}+1]
となってましたが、
αm = {-0.1(V+35)}/[exp{-(V+35)/10}-1]
αn = {-0.01(V+50)}/[exp{-(V+50)/10}+1]
の間違いでした、質問しておきながら、間違った式を紹介してしまい、まことに申し訳ありませんでした。

回答No.2

stomachman

神経細胞のモデル！いいな～。面白いことやっていますねえ。
以下、お門違いの事言ってたらゴメンしてね。

●am(t),bm(t),ah(t),bh(t),an(t),bn(t)。これらは全部v(t)を与えれば決まる関数のようです。従って、I(t)とv(t)のデータが精密に得られているならば、v(t),am(t),bm(t),ah(t),bh(t),an(t),bn(t)は全て既知ということになり、
dm(t)/dt=am(t)-(am(t)+bm(t)) m(t)
dh(t)/dt=ah(t)-(ah(t)+bh(t)) h(t)
dn(t)/dt=an(t)-(an(t)+bn(t)) n(t)
これらはm(0),h(0),n(0)が分かれば数値的に解ける。

●仮に解けたとしますとm(t),h(t),n(t)が与えられ、I(t)とv(t)も分かっているんだから、
C dv(t)/dt=I(t)-Gnamax (m(t)^3)h(t)(v(t)-Ena)-Gkmax (n(t)^4)(v(t)-Ek)-Gl (v(t)-El)
これは
C dv(t)/dt=I(t)-Gnamax (m(t)^3)h(t)v(t)+Gnamax Ena (m(t)^3)h(t)-Gkmax (n(t)^4)v(t)+Gkmax Ek (n(t)^4)-Gl v(t)+Gl El
つまりC, Gnamax, (Gnamax Ena),Gkmax,(Gkmax Ek),Gl ,(Gl El)という８個の定数を推定する問題になります。dv(t)/dtになるとさすがにノイズの影響が出てくる、ということなら、両辺をそれぞれ項別に積分しちゃえば良い。（v(0)=0は分かっています。）
ε(t) = -C v(t)+intI(t)-Gnamax int[(m(t)^3)h(t)v(t)]+Gnamax Ena int[(m(t)^3)h(t)]-Gkmax int[(n(t)^4)v(t)]+Gkmax Ek int[(n(t)^4)]-Gl int[v(t)]+Gl El t
これってε(t) の適当な評価（二乗積分とか）を最小化する線形のfitting問題です。

●従って再び
dm(t)/dt=am(t)-(am(t)+bm(t)) m(t)
これを解けるか？という問題。
m(0)を与えれば計算出来そうですから、m(0), h(0), n(0)という３個のパラメータを探索して、ε(t) が小さくできるやつを探すというのでどうでしょうか。

回答No.1

stomachman

知りません。知りませんが...
　QPchanさん自身が2つの解の候補の優劣を毎回自分で判断するとしたら、どうなさいますか？その評価方法を目的関数にすれば良い訳です。閻魔様になったつもりで考えるんですね。
　どうやら求めるのはx-y平面に描けるグラフで、パルスの波形(yの誤差)とパルスのタイミング(xの誤差)が問題のようですが、x,yはそれぞれ自由にスケールを設定できるため、それぞれにどれだけの精度を要求するかを指定しなくてはどうにもなりません。両方がそれぞれ一定の精度を満たすような解を求めたい訳ですが、GAを使うなら、その精度に満たないものの優劣を判定できなくてはなりません。
　すなわち
(1)「重み付けして和を取るっていう方法」
(2)「一方を制約条件として考える」
以外に本質的に良い方法はないと思われます。実際には(2)を適用する場合でも、条件を陽に満たす解が容易に作れない時には、制約条件をペナルティ関数で置き換えて(1)に帰着させる方法が使われ、無理に(2)をやって非線形性を強くするより良い結果が得られることが多いです。
　ただし、改良の過程で重みを変化させていくことは可能で、極端な例としてはgenerationごとに交互に一方の条件の重みを０にしてしまう、ということもできるでしょう。

　GAにこだわる理由がないのなら、もし宜しかったら具体的に微分方程式を補足してみて戴けないでしょうか。手も足も出ないおそれは多分にありますが、ひょっとしたらスマートに解決できる可能性もあります。どんな非線形性を持つ微分方程式なのか、探索空間の次元（決めたい係数の数）は幾らか、に依りますが、非線形性が強くなくて次元が低ければ、或いは。

補足

stomachmanさん、レスありがとうございます。
実際に僕がパラメータ推定をしたい微分方程式は、Hodgkin-Huxleyモデルという神経細胞のモデルです。
http://ekt715.yz.yamagata-u.ac.jp/
↓
北嶋研メモ
↓
Hodgkin-Huxley
のページに実際の微分方程式が紹介されています。

C*dv/dt=I-gnamax*m^3*h*(v-Ena)-gkmax*n^4*(v-Ek)-gl*(v-El)
dm/dt=am(v)*(1-m)-bm(v)*m
dh/dt=ah(v)*(1-h)-bh(v)*h
dn/dt=an(v)*(1-n)-bn(v)*n
am(v),bm(v),ah(v),bh(v),an(v),bn(v)は上述のページの通りです。
gnamax,gkmax,gl,Ena,Ek,El,Cは定数
v,m,h,nが状態変数
入力がI,出力がv

まず実験データ(実際に測定可能なのはvだけです)を用いてパラメータを推定する前に、シミュレーションで計算したvを用いて、gnamax,gkmax,glを最急降下法やsimulated annealingで推定しました(３個のパラメ－タを推定)、これは可能でした。

nの定常状態（ 以下n_inf(v) ）は、しばしばsigmoid関数で代用されることがあります。
n_inf(v)=1/[1+exp(-(v-n_half)/n_slope)]・・・＊
そこで、n_inf(v)を＊で表すと、
dn/dt=(n_inf-n)/tau_n・・・tau_nは時定数
dn/dt=an(v)*(1-n)-bn(v)*nをdn/dt=(n_inf-n)/tau_n代用しgnamax,gkmax,gl
,n_half,n_slope,tau_n,n^xのxを未知のパラメータとして推定してみると(７個のパラメータを推定)、最急降下法やアニ－リングなどでは、うまく行かずGAを用いることにしました。GAで推定可能でした。（最初に、ある条件で個体を絞り込み、その後パルス波形の誤差(yの誤差)を目的関数に利用。）

さらに、推定するパラメータを増やすと(dm/dt,dh/dtをsigmoidで代用して、gnamax,gkmax,gl,m_half,m_slope,tau_m,m^x1x1,h_half,h_slope,tau_h,h^x2のx2を未知のパラメータとして推定してみると(１１個のパラメータを推定)、これが７個の推定の場合と同じ目的関数を使うと、間違った値を推定します。
実際のデータに適用しても、ある程度の数のパラメータの推定までは、妥当と考えられる推定結果が得られるのですが、推定する数が多くなるとダメです。
この推定が可能ならイオンチャンネル１つを丸ごと推定できることになり、実験データの解析にも有用だと考えています。

そこで、パルスの波形(yの誤差)とパルスのタイミング(xの誤差)という２つの指標を用いる方法を探しています。一学生でパラメータ推定についも、それほど長く研究しているわけではないので、多々甘い部分はあると思います(笑)。

非線形力学系の話も全然詳しくなく、もしこの問題に非線形力学系からのアプローチで、抜け道があるのなら、ぜひ御教授していただきたいと思います。