時間tに対する1chデータ列yがありまして、それを
y=a exp(b t) + c
に対して客観的に、できれば自動的にフィッティングして、a,b,cを求めたいです。
これがただの1次関数の最小二乗法ならわかりますし、cが既知なら1次関数の応用で、というところまでもわかります。恥ずかしながら渡井には非線形最小二乗法を一般論で理解して解けるような気がしません。
Excelを使った最小二乗法手順説明サイト
http://szksrv.isc.chubu.ac.jp/lms/lms2.html
のような方法か、
C/C++のプログラム
http://www.sist.ac.jp/~suganuma/kougi/other_lect …
のようなアルゴリズムの説明をいただけると大変ありがたいです。
よろしくお願いします。
No.1ベストアンサー
- 回答日時:
回答が寄せられていませんので考えてみました。
ヒントになれば幸いです。
最小二乗法は直線に当てはめる場合の方法です。
y-c=aexp(bt)
の両辺の対数を取ると
ln(yーc)=ln(a)+bt
です。cが分かればできるのだがと書かれているのはこの式のことですね。
ならばcを推定する手順を入れたらと考えました。
時系列t1、t2、t3、・・・が等間隔とします。
ln(y2-c)-ln(y1-c)=b(t2-t1)=bτ
ln(y3-c)-ln(y2-c)=b(t3-t2)=bτ
ln(y4-c)-ln(y3-c)=b(t4-t3)=bτ
(y2-c)/(y1-c)=k
(y3-c)/(y2-c)=k
(y4-c)/(y3-c)=k
(yn-c)は等比数列になっているはずです。
y1、y2、y3に対して
(y2-c)/(y1-c)=(y3-c)/(y2-c)
より
c=(y1y3-y2^2)/(y1+y3-2y2)
が決まります。
y2、y3、y4に対しては
c=(y2y4-y3^2)/(y1+y4-2y3)
が決まります。
データがばらついていますからこのcもばらつきます。
さしあたり平均で推定していいのではないでしょうか。この値を計算して平均を求めるアルゴリズムはやさしいはずです。
全データでなくて一部でやっても構わないはずです。時間が等間隔でなければ等間隔の部分を抜き出してやればいいです。cの分布の幅も押さえておくといいでしょう。後でcを修正するときに必要になります。
cの値が推定できればln(y-c)とtのグラフを作って最小二乗法でa、bが決まります。もし適合の指標のような値も同時に得られるのでしたら少しcの値を動かして比べてみるといいと思います。cの分布の幅が分かっているとcの値を動かす目安になります。平均でやったのとあまり違わなければ平均でもいいということになります。
ご質問を見て考えたものです。素人っぽい考えかたです。
なるほど。ありがとうございます。
たしかにa,b抜きにcを求められそうですね。
落ち着いて考えてみると、今のところ収束するまでデータをとっているので、収束後の平均値でcを求めて(y-c)のデータ列を作り直してExcelの近似曲線機能でaやbを求めるという技を思いつきました。
やってみると(y-c)が0や負になると指数近似が選択できなくなるので、データ範囲に注意しながら設定するとパラメータが全部求まるようです。
定常作業なのでいつまでも続けるにはかったるいですが、しばらくは何とかなりそうです。
お教えいただいた方法をベースにデータが収束しなくてもどうにかなる自動算出方法を書いてみようと思います。
ありがとうございます。
No.4
- 回答日時:
最小自乗法を使うときにはいくつか留意点があるかと思います。
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
- 回答日時:
>今のところ収束するまでデータをとっているので、収束後の平均値でcを求めて(y-c)のデータ列を作り直して・・・
収束が起こる現象を見ているだろうというのに気がつきませんでした。
b<0になっているということですね。減衰曲線です。
でも考えてみれば当然ですね。b>0の発散の場合、測定範囲が広すぎて対応できないはずです。
式は初めから
y=aexp(-bt)+c
として考えるほうがいいですね。
ふと思ったのですが指数関数的に収束するという事に対する裏付けはあるのでしょうか。これは着目している現象の中での測定量の性質として決まってくるはずのものですね。一応確めておく必要があると思います。
収束する関数は指数関数だけでは有りません。データのばらつきが大きければどういうカーブにでも合わせることが出来ます。どういう式に従って収束するかは理論的に予測しておく必要があります。
1/t^(n)で減衰する場合は
y=a/t^(n)+c
です。
log(y-c)=loga-nlogt
となります。
手作業でやるときは半対数方眼紙を使うか両対数方眼紙を使うかの違いになります。いきなり機械に計算させてしまうのでははなくて一度図示するという判断のプロセスを入れるというのも大事なことでしょう。「自動計算」を考えておられますのでちょっと気になりました。
本当に指数関数的なのか、ってところの確認でGWを使い切っちゃたのですが、どうやら極初期以外は指数関数的って理論がたちまして、実験結果からも合理的なパラメータをはじき出すことができました。
結局大量にExcelのフィッティング機能連発という人海戦術で、グラフ書きながら確認もしました。
いろいろありがとうございます。
No.2
- 回答日時:
#1です。
途中ミスタイプがあるのに気がつきました。
念のため、訂正しておきます。
(誤)「y2、y3、y4に対しては
c=(y2y4-y3^2)/(y1+y4-2y3)
が決まります。」
(正)「y2、y3、y4に対しては
c=(y2y4-y3^2)/(y2+y4-2y3)
が決まります。」
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
似たような質問が見つかりました
- 数学 数学の問題の解き方を教えてください! 3 2022/11/02 17:32
- JavaScript 最小二乗法 2 2023/01/01 20:57
- 数学 すべての自然数とすべての実数を1対1で対応させる(すべての実数を一列に並べる)方法について 3 2023/05/26 17:14
- 数学 8 件の住宅について, 駅からの徒歩時間 (分) と賃料 (万円) を調べたところ, (徒歩時間, 1 2022/12/18 18:09
- 数学 8 件の住宅について, 駅からの徒歩時間 (分) と賃料 (万円) を調べたところ, (徒歩時間, 2 2022/12/18 20:26
- 数学 『弧は弦より長し』 8 2022/04/18 10:23
- 数学 回答の意味について 3 2023/07/06 14:14
- 数学 【高1 数学Ⅰ 二次関数】 二次関数 f(x)=x^2-4ax+8a がある。ただし、aは正の定数と 3 2022/07/23 15:46
- 数学 ラグランジュの未定乗数法を用いる問題 3 2023/05/15 14:48
- 数学 数学?算数の問題です どのような解答になりますか? 2 2022/04/22 04:46
このQ&Aを見た人はこんなQ&Aも見ています
おすすめ情報
このQ&Aを見た人がよく見るQ&A
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
gnuplot でこのような濃淡グラ...
-
FFTのDC成分って、なんで大きく...
-
隣接平均と移動平均
-
指数関数のカーブフィッティング
-
近似直線について
-
MOSFETのgm-Vgs特性
-
極点図のプロット用ソフトウェア
-
B-H曲線について
-
MTFを求める際に、平滑化や隣接...
-
フーリエ級数展開を使ってどう...
-
1/3オクターブバンドについて質...
-
静的・動的の意味
-
Igorの使い方について
-
どなたか、S45Cの電気伝導率(S...
-
10,000百万円っていくらですか?
-
時間を100進法であらわしたい。
-
「強度」は高い?強い?
-
合成関数の微分を使う時と、使...
-
穴あけのあけってどんな漢字で...
-
「1人あたりの1年間あたりの~...
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
おすすめ情報