最新閲覧日:

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

A 回答 (2件)

データを(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のひとつの山の一部だけのデータしかなければ、どの山のものなのか特定するのは極めて困難である。
●計算を1回やれば良いのか、数十回やるのか、どんどん入ってくるデータを自動処理したいのか。
●a,b,cの大体の値は分かっているのか。
などです。これらの条件を勘案して、適当な方法を決めなくてはならない。ともかく幾つか汎用性の高いアプローチを示しましょう。
教科書に載っている「非線形最小二乗法」は、まあ行ってみれば最後の手段です。
(でも「最小二乗法による実験データの解析」東大出版会 は持っていて損のない本です。)

★アプローチ(1)
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を計算することを繰り返せばよい。

★アプローチ2
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次、1次を計算すれば解けて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
これを強引に項別に2回部分積分すると
(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を非線形最小二乗法で改良すると、出発値が良いので非常に安定かつ短時間で計算できます。(ガウス-ニュートン法程度で十分。)

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

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

お礼日時:2001/01/29 00:39

x,y が測定値で,a,b,c を最小2乗法で決めたいんですね.


で,a,b,c について線型でないから困っている,ですね.

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

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

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

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

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

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

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

お礼日時:2001/01/29 00:36

お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!

このQ&Aを見た人はこんなQ&Aも見ています

このQ&Aを見た人が検索しているワード


このQ&Aを見た人がよく見るQ&A

人気Q&Aランキング

おすすめ情報