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を非線形最小二乗法で改良すると、出発値が良いので非常に安定かつ短時間で計算できます。（ガウス-ニュートン法程度で十分。）

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

お礼

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

回答No.1

siegmund

x,y が測定値で，a,b,c を最小２乗法で決めたいんですね．
で，a,b,c について線型でないから困っている，ですね．

よく使われる方法は２つ
いずれも，あらっぽくて良いですから，a,b,c の値の
大体の見当をつけておきます．

(a) a,b,c を適当なステップで変えて残差の２乗の和を計算する．
残差の２乗の和が最小になる a,b,c が求めるものです．
例えば，a,b,c をどれも１０通り変えて，1000 通りについて残差の２乗和
を見ればよい．
１回ではステップが荒すぎるなら，残差の２乗和が最小になっている
あたりをもっと細かくやる．
例えば，範囲を１回目の 1/10 に絞り，ステップも 1/10，
a,b,c の値はやはり１０通りずつですね．
お好きなだけ手続きを繰り返してください．

(b) 線型じゃなくて困るなら，線型にする．
a の見当値を a0 とし，a = a0 + δa とする．
b,c についても全く同様．
で，δa, δb, δc の１次まで展開すれば，
δa, δb, δc について線型になりますね．
線型の最小２乗法だから簡単にできる．
できたら，a0 + δa を新たに見当値 a0 と思って，
同じことを繰り返す．
何回かやれば十分収束します．

定めるべき係数がべらぼうにたくさんあるときは，
収束を早くする線型代数的テクニックが使われますが，
今は３個しかありませんから，
テクニック学んだりプログラムが複雑になったりするよりは
素直にやるのが早いでしょう．

どちらも簡単なプログラムでしょう．
主要部分は繰り返しですから，
適当に判断してループから抜ければいいですね．

非常に運が悪いとうまくいかないことがまれにありますが
たいていは上のどちらかで十分です．
誤差評価も考えてくださいね．
誤差分より細かいステップまでやっても意味がありません．

お礼

回答ありがとうございます。
参考にさせていただきます。