指数関数のカーブフィッティング

質問

時間ｔに対する１ｃｈデータ列ｙがありまして、それを
ｙ＝ａ　ｅｘｐ（ｂ　ｔ）　＋　ｃ
に対して客観的に、できれば自動的にフィッティングして、ａ，ｂ，ｃを求めたいです。
これがただの1次関数の最小二乗法ならわかりますし、ｃが既知なら１次関数の応用で、というところまでもわかります。恥ずかしながら渡井には非線形最小二乗法を一般論で理解して解けるような気がしません。

Excelを使った最小二乗法手順説明サイト
http://szksrv.isc.chubu.ac.jp/lms/lms2.html
のような方法か、
C/C++のプログラム
http://www.sist.ac.jp/~suganuma/kougi/other_lecture/SE/predict/predict.htm#2
のようなアルゴリズムの説明をいただけると大変ありがたいです。

よろしくお願いします。

回答No.1

htms42 · ベストアンサー

回答が寄せられていませんので考えてみました。
ヒントになれば幸いです。

最小二乗法は直線に当てはめる場合の方法です。
ｙ－ｃ＝ａｅｘｐ（ｂｔ）
の両辺の対数を取ると
ｌｎ（ｙーｃ）＝ｌｎ（ａ）＋ｂｔ
です。ｃが分かればできるのだがと書かれているのはこの式のことですね。
ならばｃを推定する手順を入れたらと考えました。

時系列ｔ１、ｔ２、ｔ３、・・・が等間隔とします。
ｌｎ（ｙ２－ｃ）－ｌｎ（ｙ１－ｃ）＝ｂ（ｔ２－ｔ１）＝ｂτ
ｌｎ（ｙ３－ｃ）－ｌｎ（ｙ２－ｃ）＝ｂ（ｔ３－ｔ２）＝ｂτ
ｌｎ（ｙ４－ｃ）－ｌｎ（ｙ３－ｃ）＝ｂ（ｔ４－ｔ３）＝ｂτ

（ｙ２－ｃ）／（ｙ１－ｃ）＝ｋ
（ｙ３－ｃ）／（ｙ２－ｃ）＝ｋ
（ｙ４－ｃ）／（ｙ３－ｃ）＝ｋ

（ｙｎ－ｃ）は等比数列になっているはずです。
ｙ１、ｙ２、ｙ３に対して
（ｙ２－ｃ）／（ｙ１－ｃ）＝（ｙ３－ｃ）／（ｙ２－ｃ）
より
ｃ＝(ｙ１ｙ３-ｙ２＾２）／（ｙ１＋ｙ３－２ｙ２）
が決まります。
ｙ２、ｙ３、ｙ４に対しては
ｃ＝(ｙ２ｙ４-ｙ３＾２）／（ｙ１＋ｙ４－２ｙ３）
が決まります。
データがばらついていますからこのｃもばらつきます。
さしあたり平均で推定していいのではないでしょうか。この値を計算して平均を求めるアルゴリズムはやさしいはずです。
全データでなくて一部でやっても構わないはずです。時間が等間隔でなければ等間隔の部分を抜き出してやればいいです。ｃの分布の幅も押さえておくといいでしょう。後でｃを修正するときに必要になります。

ｃの値が推定できればｌｎ（ｙ－ｃ）とｔのグラフを作って最小二乗法でａ、ｂが決まります。もし適合の指標のような値も同時に得られるのでしたら少しｃの値を動かして比べてみるといいと思います。ｃの分布の幅が分かっているとｃの値を動かす目安になります。平均でやったのとあまり違わなければ平均でもいいということになります。

ご質問を見て考えたものです。素人っぽい考えかたです。

お礼

なるほど。ありがとうございます。
たしかにａ，ｂ抜きにｃを求められそうですね。

落ち着いて考えてみると、今のところ収束するまでデータをとっているので、収束後の平均値でｃを求めて（ｙ－ｃ）のデータ列を作り直してExcelの近似曲線機能でａやｂを求めるという技を思いつきました。

やってみると（ｙ－ｃ）が０や負になると指数近似が選択できなくなるので、データ範囲に注意しながら設定するとパラメータが全部求まるようです。
定常作業なのでいつまでも続けるにはかったるいですが、しばらくは何とかなりそうです。
お教えいただいた方法をベースにデータが収束しなくてもどうにかなる自動算出方法を書いてみようと思います。

ありがとうございます。

回答No.4

foobar

最小自乗法を使うときにはいくつか留意点があるかと思います。

y=Aexp(Bx)に最小自乗法でフィッティングさせるとき、
y=Aexp(Bx)をそのまま使う場合と、
log y=logA+Bxの形にして最小自乗法を適用するときで、
AとBの値に違いができます。
logyで計算すると、log y とフィッティング曲線との差が均等化されるため、yが大きなところでの誤差が（実数のグラフで見ると）大きくなります。
特に、減衰してゆくデータで、十分減衰して0付近になったデータが多いと、A,Bの推定値がここの部分に引っ張られてしまう可能性があります。

どういうデータに適用するか、推定したA,Bをどのように使うかにも寄りますが、場合によっては、直接 y=Aexp(Bx)の形で、数値計算的にA,Bを求めるほうがよい場合があります。

たとえば、
a)A,B,Cの値の初期値を決める（Cは十分減衰したところのデータの平均値、A,Bはlog(y-C)=logA+Bxの形で最小自乗法により決定）
b)e0=Σ(y0-y)^2を計算　（y0はAexp(Bx)+Cの値）
c)Aを微小量(正負)変えて、同様にe1,e2を計算。
d)e0,e1,e2から誤差が最小になるAの値を推定
e)B,Cについても同様に新しい値を推定
f)A,B,Cが所定の誤差内に収まるまで、b)からe)を繰り返す。
みたいな手法もあります。（最初の初期値のとり方がまずいと変な値を出す可能性があるとか、dの推定の手順がまずいと数値が振動して収束しない、計算終了の判断条件をどうするか、といった問題点がありますが）

お礼

今回の測定で誤差要因が何なのかが正しく評価できていないと、対数とってから最小自乗法するのがよいのか、もとの値のまま誤差を最小にするのがよいのか判断がつきません。そして、今のところそこをまじめに考えていくほどの余裕もないのでテキトーに誤魔化したい気持ちでいっぱい。実際、今回使ったExcelが、どういうアルゴリズムで指数フィッティングを実施しているのかも不明。参考値にはなっても、まじめな値にはならないのを覚悟で、今年度後半のまじめな細かい研究の指針をオーダー計算する補強実験としては、なんとかなった模様。

いつまでもcsvで書き出して、Excelでせっせと計算して、グラフ書いてフィッティングしてパラメータ表示させて、手で書き留める、をいくつも繰り返すなんてことをやってられませんので、夏くらいまでにはご指摘のアルゴリズムも候補の一つとして検討いたします。

ご教授いただきありがとうございました。

回答No.3

htms42

＞今のところ収束するまでデータをとっているので、収束後の平均値でｃを求めて（ｙ－ｃ）のデータ列を作り直して・・・

収束が起こる現象を見ているだろうというのに気がつきませんでした。
ｂ＜０になっているということですね。減衰曲線です。
でも考えてみれば当然ですね。ｂ＞０の発散の場合、測定範囲が広すぎて対応できないはずです。
式は初めから
ｙ＝ａｅｘｐ（－ｂｔ）＋ｃ　　
として考えるほうがいいですね。

ふと思ったのですが指数関数的に収束するという事に対する裏付けはあるのでしょうか。これは着目している現象の中での測定量の性質として決まってくるはずのものですね。一応確めておく必要があると思います。
収束する関数は指数関数だけでは有りません。データのばらつきが大きければどういうカーブにでも合わせることが出来ます。どういう式に従って収束するかは理論的に予測しておく必要があります。

１／ｔ^(n)で減衰する場合は
ｙ＝ａ／ｔ^(n)＋ｃ　
です。
ｌｏｇ（ｙ－ｃ）＝ｌｏｇａ－ｎｌｏｇｔ
となります。

手作業でやるときは半対数方眼紙を使うか両対数方眼紙を使うかの違いになります。いきなり機械に計算させてしまうのでははなくて一度図示するという判断のプロセスを入れるというのも大事なことでしょう。「自動計算」を考えておられますのでちょっと気になりました。

お礼

本当に指数関数的なのか、ってところの確認でＧＷを使い切っちゃたのですが、どうやら極初期以外は指数関数的って理論がたちまして、実験結果からも合理的なパラメータをはじき出すことができました。
結局大量にExcelのフィッティング機能連発という人海戦術で、グラフ書きながら確認もしました。

いろいろありがとうございます。

回答No.2

htms42

＃１です。

途中ミスタイプがあるのに気がつきました。
念のため、訂正しておきます。

（誤）「ｙ２、ｙ３、ｙ４に対しては
ｃ＝(ｙ２ｙ４-ｙ３＾２）／（ｙ１＋ｙ４－２ｙ３）
が決まります。」

（正）「ｙ２、ｙ３、ｙ４に対しては
ｃ＝(ｙ２ｙ４-ｙ３＾２）／（ｙ２＋ｙ４－２ｙ３）
が決まります。」