ガウスフィッティングについて

時間 t に対するデータ列 yi (i=1, 2, …) に
最小二乗法を用いて
ガウス関数 y(t) = a + b exp( -(t-c)^2 / 2σ^2 ) を
当てはめたいのです。
これまでは、エクセルなどのソフトを用いて
未知数 a, b, c, σを求めていたのですが、
これをソフトなどを用いずに理論式で解くにはどうすれば良いでしょうか？

ひとつのやり方として、
ガウス関数の右辺の a を左辺に移行して、両辺の対数を取り
対数データ列 ln(yi - a) との誤差を求めて最小二乗法…
という手を考えたのですが、
これだと、誤差の総和を求める際、Σの中にパラメータの a が
ln(yi - a) という形で残ってしまい、
誤差の偏微分で得られるパラメータに関する方程式が、
非線形な連立方程式となって解くことができません。
(説明が分かり難くて申し訳ありません)

その他にもいろいろ考えたのですが、
どうにもうまくいきません。
とにかく、パラメータの a が邪魔で、
これを b, c, σを用いずに先に求める、
という方向で考えているのですが、
何か良い方法がありましたら、よろしくお願いします。

No.1 ベストアンサー 回答者： stomachman 回答日時： 2011/03/02 15:09

　最小化したいのは

E=Σ(yi-y(t))^2
ですね。パラメータをベクトルで
p=(a,b,c,σ)
と書くことにします。

　どうやっても非線形で、反復計算が必要です。「理論式」という観点で言い換えるなら、「最小二乗解pに収束する漸化式」を作るところまでしか出来ない、ってことです。
　その漸化式は、たとえば「第n回目の反復で得られたパラメータの推定値p[n]を使って、(∂E/∂p)(p[n])を係数とするpの一次式でEを近似し、これを使ってp[n]よりも最小二乗解に近いp[n+1]を計算する」というものですから、「適切な初項を与えれば最小二乗解に収束する漸化式」でしかありません。（他にもアプローチはありますが、初項が必要なのは同じ。）
　ここで「適切な初項」というのは、「初項で決まるガウス曲線がもしデータとかけ離れていたら、漸化式が収束しない」ってことです。では初項が「適切」かどうかを予め判定する手段はというと、これまた一般にはありませんで、計算してみて収束していくかどうかを見るしかない。そういう事情です。

　なので、実地に計算するには、漸化式の初項（つまりpの最初の推測値）を決める方法が必要。これはイロイロ工夫するしかありません。

　以下は「適切な初項を決める手段」の話なので、現段階ではご興味ないのかも知れませんが：

　「パラメータaをなんとか推定したうえで、対数を取って線形最小二乗法でb,c,σを出す」というのは一つの方法です。ただしこのままだと
F=Σ(ln(yi)-ln(y(t)))^2
を最小化することになり、Eの最小化とはかなり違います。むしろ
G=Σ(yi (ln(yi)-ln(y(t)))^2
でEをよく近似できます。( ln(1+x)≒x を使っています。）
　もちろん、これは飽くまで初項を推定する方法であって、それを使って「Eを最小化する漸化式」で反復計算をすることになります。

　別の手段として、もしデータyiによくフィットする曲線y(t)が存在する（データがほぼ綺麗にガウスカーブになっている）なら、データの最大値付近に放物線をフィッティングするなどしてcの推測値を決め、さらに直線
y=(最大値-最小値)/2
とデータの折れ線グラフとの交点の間隔2w（半値幅）からσを推定して、線形最小二乗法でa,bを決める、などなど。

この回答への補足

> 「最小二乗解pに収束する漸化式」を作るところまでしか出来ない
いろいろ調べているのですが、やはり反復法などが多くヒットしますね。
やはり、理論的に解くというのは厳しいものなのか…。
その「漸化式を作る」というのを、もう少し詳しく解説して貰えないでしょうか？

それから、E を最小化する p と F を最小化する p は違う解になるのですか？
よく理解していませんので、自分的には同じになると思っておりましたが。

補足日時：2011/03/03 12:50

この回答へのお礼

ご回答どうもありがとうございました。
頂いた回答を参考にしつつ、いろいろ調べた結果、
なんとか、解法に辿り着きました。
どうもありがとうございました。

お礼日時：2011/03/07 18:48

No.2 回答者： spring135 回答日時： 2011/03/02 15:22

必ずしも良く理解しているわけではないのですが確率密度関数としてガウス分布を仮定するならば

a=0, c=Σyi/n, b=1/(2π)^.5σ、σ=Σ(yi-c)^2/n
でデータから一義的に来ますと思いますが、データセットが違うのでしょうか。

この回答への補足

データは、その両端において必ずしも 0 に収束しません。
両端では、ある値(正か負)付近で振動しているようなデータです。
ですので、a は必ずしも 0 になりません。

補足日時：2011/03/03 12:22

この回答へのお礼

ご回答どうもありがとうございました。
頂いた回答を参考にしつつ、いろいろ調べた結果、
なんとか、解法に辿り着きました。
どうもありがとうございました。

お礼日時：2011/03/07 18:48