データ列の曲線によるフィッティング

データ列(xi,yi)(i=1,...,n)を関数
y=α(x+γ)^2 exp(-β(x+γ)) (α>0,β>0,γ>0)
でフィッティングしたいです。対数をとって普通に最小２乗法で解こうとしたら得られる連立方程式が線形でなくて解けませんでした。
どうしたらいいのでしょうか？

No.1 ベストアンサー 回答者： stomachman 回答日時： 2002/01/15 11:33

完全な意味での最小二乗法ではないけれど、無理矢理線形の問題に帰着して簡単に計算するstomachman我流のやりかたがあります。

（超平面法と勝手に呼んでいます。）この方法は、「データがモデルに良く合う」という時に限り「最小二乗解とほぼ一致する」という性質があります。

モデルをデータと区別して、
f(t) = α((t+γ)^2) exp(-β(t+γ))
と書きましょう。
両辺をtで微分すると
f'(t) = 2α(t+γ) exp(-β(t+γ))-βα((t+γ)^2) exp(-β(t+γ))
(t+γ)f'(t) = 2α(t+γ)^2) exp(-β(t+γ))-(t+γ)βα((t+γ)^2) exp(-β(t+γ))
(t+γ)f'(t) = 2f(t)-(t+γ)βf(t)
(t+γ)f'(t) = (2-βγ)f(t)-βtf(t)
という微分方程式が得られます。
　この両辺をtで積分すると、(t=x1～xiまで）
∫(t+γ)f'(t)dt = (2-βγ)∫f(t)dt-β∫tf(t)dt
左辺を部分積分して
(xi+γ)f(xi)-(x1+γ)f(x1)-∫f(t)dt = (2-βγ)∫f(t)dt-β∫tf(t)dt
整理すると
xif(xi)-x1f(x1) = (3-βγ)∫f(t)dt-β∫tf(t)dt-γ(f(xi)-f(x1))
です。
ここで、f(xi)の代わりにyiを使っても、「データとモデルが良く合う」のなら大差は出ない筈。だから
∫f(t)dtの代わりにyをx1～xiまで数値積分したものJi、
∫tf(t)dtの代わりにxyをx1～xiまで数値積分したものKi、を用いることにして、残差をεiと書くことにし、
xiyi-x1y1 = (3-βγ)Ji-βKi-γ(yi-y1)+εi　(i=1,2,....,n)

改めて
Yi = xiyi-x1y1
Li = yi - y1
A = (3-βγ)
B = -β
C = -γ
とおくと、
Yi = A Ji + B Ki + C Li + εi　(i=1,2,....,n)
という線形問題になったわけで、これは簡単に解けます。
　ところで、A,B,Cはβとγだけで表される一方、αが出てきませんね。
つまり
・A,B,Cからβ,γを求めると１通りに決まりません。これは気にしないで、たとえばAを無視してβ,γを計算します。ついでに、こうして得たβ,γを使って(3-βγ)を計算してみます。もしこれがAとまるで合わないようなら、「データがモデルに良く合う」という条件が満たされていない。yiに含まれるノイズが大きすぎるんです。
・αを求めるには、
f(t) = α((t+γ)^2) exp(-β(t+γ))
を再び使います。γとβは既に分かっているから、簡単な線形最小二乗法の問題ですね。

より高精度の最小二乗解を求めるには、超平面法で得た解を初期値にして、非線形最小二乗法（ガウス・ニュートン法など）を使って改良すると良いです。（教科書としてはいつも、中川・小柳「最小二乗法による実験データ解析」東京大学出版会をお薦めしています。）

No.3 回答者： noname#21649 回答日時： 2002/01/16 07:33

最小二乗法と非線型２情報の話しは既に出ていますので．ＯＲの話しにしましょう。

最適化問題で非線型問題がどうしても融けない場合に使います。
まず．でたらめにアルファ・べーた・ガンマーの値を決めます。
簡単にするのでしたらば対数で１ｅ－８から１え＋８の間が均等になるような点（１．１０．１００とか標準数（ＪＩＳＺ番号忘れ）とかあります）を選び．残差を求めます。
そのなかで．もっとも残差が少ない．数値の組み合わせを選びます。

次に．どれか一つをちょっと変えて計算してみます。「ちょっと」の程度ですが．切りの良い（２進数で考えてください）１／１２８ぐらいではじめます。もし．最初の計算よりも残差が少なくなったらばその解をえらびます。もし．残差がおおきくなってしまったらは．「ちょっと」の程度を更に１／２にします。計算していて．計算機Ｅの問題が出てきたら計算を止めます。
この考え方を基にして．３つの変数の解の付近（目安としては前回の計算した範囲程度）を１２８分割して演算をしてその中のもっとも残差の少ない点を解とします。

最低でも倍精度演算が必須です。誤差関数を求めるために８倍精度とか無限倍精度（文字演算で行うために計算がやたら遅いのでやりたくない）とかを使う方法もあります。

ＯＲ問題の亜種です。私は場当たり法と呼んでいます。メッシュで計算する法ですから．解の巣問題や鞍部問題（１の方の参考文献を読んでください）に対して（計算を間違えなければ）強力な計算方法です。しかし．計算時間がやたらかかるので．評判はよろしくないです。８－１０変数で．１００点．８０４８６で８時間なんてざらです（真の最適解とは限らないので．初期値をちょっとずらして再計算するとすぐに１－２ヶ月使ってしまいます）から。

No.2 回答者： siegmund 回答日時： 2002/01/15 11:43

これは非線形最小二乗法の問題になりますが，

過去いくつかのスレッドがあり，処方もいろいろ回答されています．
たとえば，
http://oshiete1.goo.ne.jp/kotaeru.php3?q=179919
http://oshiete1.goo.ne.jp/kotaeru.php3?q=190632
http://oshiete1.goo.ne.jp/kotaeru.php3?q=98446
http://oshiete1.goo.ne.jp/kotaeru.php3?q=97271
http://oshiete1.goo.ne.jp/kotaeru.php3?q=33318
など．
他にも質問検索で「最小二乗法」「最小自乗法」とやるとかなりヒットします．