-PR-
締切り
済み

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

  • すぐに回答を!
  • 質問No.33318
  • 閲覧数1697
  • ありがとう数10
  • 気になる数0
  • 回答数2
  • コメント数0

お礼率 100% (2/2)

データ配列に、
y = a sinc(b(x-c))
で表せれる式をフィッティング(最小二乗法など)したいのですが、良い方法がわかりません。
どなたか教えてもらえませんでしょうか?
通報する
  • 回答数2
  • 気になる
    質問をブックマークします。
    マイページでまとめて確認できます。

回答 (全2件)

  • 回答No.1
レベル13

ベストアンサー率 64% (700/1089)

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 ...続きを読む
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個しかありませんから,
テクニック学んだりプログラムが複雑になったりするよりは
素直にやるのが早いでしょう.

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

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

お礼率 100% (2/2)

回答ありがとうございます。
参考にさせていただきます。
投稿日時 - 2001-01-29 00:36:40


  • 回答No.2
レベル14

ベストアンサー率 57% (1014/1775)

データを(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]をたとえば最小二乗法で小さく押さえ込みたい。という問題ですね。 いろんな方法があるんです。そして、アプローチを決める上でのポイントは ●s ...続きを読む
データを(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を非線形最小二乗法で改良すると、出発値が良いので非常に安定かつ短時間で計算できます。(ガウス-ニュートン法程度で十分。)

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

お礼率 100% (2/2)

丁寧な回答ありがとうございます。
早速回答中のいくつかのアプローチで考えてみます。
わからないことがあったら、また質問させていただきます。そのときはまたよろしくお願いします。
投稿日時 - 2001-01-29 00:39:25
このQ&Aで解決しましたか?
関連するQ&A
-PR-
-PR-
このQ&Aにこう思った!同じようなことあった!感想や体験を書こう
このQ&Aにはまだコメントがありません。
あなたの思ったこと、知っていることをここにコメントしてみましょう。

その他の関連するQ&A、テーマをキーワードで探す

キーワードでQ&A、テーマを検索する
-PR-
-PR-
-PR-

特集


いま みんなが気になるQ&A

関連するQ&A

-PR-

ピックアップ

-PR-
ページ先頭へ