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

質問

データ配列に、
y = a sinc(b(x-c))
で表せれる式をフィッティング（最小二乗法など）したいのですが、良い方法がわかりません。
どなたか教えてもらえませんでしょうか？

回答No.2

stomachman

データを(x[i],y[i]) (i=1,2,....N)とし、x[i]は誤差が無視でき、y[i]には誤差がある。つまり、
f(a,b,c;x) = a sinc(b(x-c))
としたとき、
y[i]-f(a,b,c:x) = ε[i]
このε[i]をたとえば最小二乗法で小さく押さえ込みたい。という問題ですね。

いろんな方法があるんです。そして、アプローチを決める上でのポイントは
●sincが波打ってる所まで、つまり遠くまでxがあるのかどうか。
●xが沢山あるかどうか。波の一つあたり少なくとも幾つかはxがあるのかどうか。
●xが等間隔に与えられているのか。
●yのノイズがモデルに比べてものすごく大きかったりしないか。
●xにノイズはないか。
●sincのひとつの山の一部だけのデータしかなければ、どの山のものなのか特定するのは極めて困難である。
●計算を１回やれば良いのか、数十回やるのか、どんどん入ってくるデータを自動処理したいのか。
●a,b,cの大体の値は分かっているのか。
などです。これらの条件を勘案して、適当な方法を決めなくてはならない。ともかく幾つか汎用性の高いアプローチを示しましょう。
教科書に載っている「非線形最小二乗法」は、まあ行ってみれば最後の手段です。
（でも「最小二乗法による実験データの解析」東大出版会　は持っていて損のない本です。）

★アプローチ（１）
sin(A+B) = sinAcosB+cosAsinB
を使って、
f(a,b,c;x) x=c f(a,b,c;x)+{(a/b)(cos bc)} (sin bx) +{-(a/b) (sin bc)} (cos bx)
とモデルを書き換えます。
f(a,b,c;x)=y[i]-ε[i]
を使って、
y[i](x[i])=
　　c y[i]+
　　{(a/b)(cos bc)} (sin bx[i]) +
　　{-(a/b) (sin bc)} (cos bx[i])+(x[i]-c)ε[i]
ですね。従って、
Y[i]= y[i](x[i])
α=c
β={(a/b)(cos bc)}
γ={-(a/b) (sin bc)
δ[i]=(x-c)ε[i]
と考えてしまえば、これはbだけが非線形で、あとのa,cについては線形最小二乗法の問題です。だからbを適当な方法で探索すれば良い。その過程で一旦cの概算値が分かってしまえば、
y[i](x[i])/(x[i]-c)=
　　c y[i]/(x[i]-c)+
　　{(a/b)(cos bc)} (sin bx[i])/(x[i]-c) +
　　{-(a/b) (sin bc)} (cos bx[i])/(x[i]-c)+ε[i]
によって、より良い評価が得られるから、bを決めてはc,aを計算することを繰り返せばよい。

★アプローチ２
xが等間隔でたっぷりあって、yのノイズが小さいのなら、フーリエ変換してしまいましょう。
F(a,b,c;ω) = (a/√(2π) )integral{x=-∞～∞} sinc(b(x-c)) exp[-iωx] dx
とすると
F(a,b,c;ω) = (a/√(2π) )integral{x=-∞～∞} sinc(bx) exp[-iω(x+c)] dx
=exp[-iωc] (a/(b√(2π) ))integral{x=-∞～∞} {sin(bx)/x} exp[-iωx] dx
=exp[-iωc] (a/(b√(2π) ))√(2π) ---- ｜ω｜>｜b｜　さもなくば0
=exp[-iωc] (a/b)---- ｜ω｜>｜b｜　さもなくば0
=(a/b) ((cos ωc)-i(sin ωc))---- ｜ω｜>｜b｜　さもなくば0
ここで
|F(a,b,c;ω|^2 =(a/b)^2---- ｜ω｜>｜b｜　さもなくば0
これでa, bを決めてしまうことができる（かも知れません）。そうすれば、
既知のものを大文字で書けば
f(A,B,c;x) x=c f(A,B,c;x)+{(A/B)(cos Bc)} (sin Bx) -(A/B) (sin Bc) (cos Bx)
なので、
y[i]x[i]=c y[i]x[i]+(A/B)(cos Bc) (sin Bx[i]) -(A/B) (sin Bc) (cos Bx[i])+(x[i]-c)ε[i]

未知数はcだけですね。この線形方程式はBが既知なのでフーリエ級数展開の0次、１次を計算すれば解けてc, (A/B)(cos Bc), (A/B) (sin Bc)が得られます。cが3通りも求められてしまう。これを使って適当にcを選んで
y[i]x[i]/(x[i]-c)=c y[i]x[i]/(x[i]-c)+{(A/B)(cos Bc)/(x[i]-c)} (sin Bx[i]) -{(A/B) (sin Bc)/(x[i]-c)} (cos Bx[i])+ε[i]
で最小二乗法をやっても良いし、あとはいろいろやり方ありますねえ。

★3アプローチ3
ちょっと乱暴ですが、一気に線形で解きます。データの誤差が少なく、データ数が多い場合には非常に強力。
モデルf(a,b,c;x)は微分方程式
(b^2)(x-c)f-f'+(x-c)f''=0
の解である。したがって、
(b^2)xf-(cb^2)f-f'+xf''-cf''=0
これを強引に項別に２回部分積分すると
(b^2)int int xf (dx)^2-(cb^2)int int f (dx)^2-int f dx + xf-int fdx -intint f (dx)^2 - cf =Cx+D+誤差
という線形問題になる。ここでfにyを代入して数値積分し（従ってどの積分も定数になる）、線形最小二乗法で誤差の二乗和を最小化すれば、
(b^2), (cb^2), c, C,D
という係数を決定することが出来る。あとのaを決める方法は自明でしょう。

こうやって得たa,b,cを非線形最小二乗法で改良すると、出発値が良いので非常に安定かつ短時間で計算できます。（ガウス-ニュートン法程度で十分。）

ちと眠たくて、推敲が甘いですから、式の間違いなどチェックはご自分で。
意味不明の部分は、ご遠慮なく補足質問してください。

お礼

丁寧な回答ありがとうございます。
早速回答中のいくつかのアプローチで考えてみます。
わからないことがあったら、また質問させていただきます。そのときはまたよろしくお願いします。