curve_fitが正弦波でも失敗する

Question

もっと複雑な問題に移る前に、私はcurve_fitを使って単純な正弦波（雑音でもない）をテストとして試そうとしています。残念ながら、それは遠隔で正しい答えを与えることさえありません。ここに私の構文です：

x = linspace(0,100,300) 
y = sin(1.759*x) 
def mysine(x, a): 
    return sin(a*x) 

popt, pcov = curve_fit(mysine, x, y) 
popt 
array([ 0.98679056]) 

そして、私は最初の推測（たとえば1.5）しようとした場合：

popt, pcov = curve_fit(mysine, x, y, p0=1.5) 
popt 
array([ 1.49153365]) 

...正しい答えほど遠いまだありません。

私は、関数がどのくらいうまくサンプリングされているかを考えれば、フィットがうまくいかないことに驚いていると思います。

Answer 1

カーブフィッティングが必ずしも単純ではありません。 curve_fitアルゴリズムは、最小二乗カーブフィッティングに基づいており、通常、入力パラメータの初期推測が必要です。あなたが適合したい機能の種類に応じて、最初の推測は良いものでなければなりません。

初期の推測を試みたにもかかわらず、サンプリング周波数とウェーブの周波数に関係する追加の問題があると思います。詳細については、WikipediaのNyquist-Shannon sampling theoremをご覧ください。簡単に言えば、あなたの波の周波数は1.759 /（2 * pi）= 0.28で、x配列（〜0.33）のサンプリング周波数に非常に近いことがわかります。起こりうるもう一つの問題は、あなたの機能に合わせて振動をあまりにも多く持つことです。

コードを動作させるためには、ウェーブの周波数を上げる（a> 4 * 0.33）か、サンプリング周波数を上げて空間ベクトルの長さを短くすることをお勧めします。x

私は、次のコードを実行し、図示hereような結果を得た：

# -*- coding: utf-8 -*- 
import numpy as np 
import pylab as pl 
from scipy.optimize import curve_fit 

def mysine(x, a): 
    return 1. * np.sin(a * x) 

a = 1.759 # Wave frequency 
x = np.linspace(0, 10, 100) # <== This is what I changed 
y = np.sin(a * x) + 0. * np.random.normal(size=len(x)) 

# Runs curve fitting with initial guess. 
popt, pcov = curve_fit(mysine, x, y, p0=[1.5]) 

# Calculates the fitted curve 
yf = mysine(x, *popt) 

# Plots results for comparison. 
pl.ion() 
pl.close('all') 
fig = pl.figure() 
ax = fig.add_subplot(111) 
ax.plot(x, y, '-', c=[0.5, 0.5, 0.5]) 
ax.plot(x, yf, 'k-', linewidth=2.0) 
ax.text(0.97, 0.97, ur'a=%.4f, ã=%.4f' % (a, popt[0]), ha='right', va='top', 
    fontsize=14, transform=ax.transAxes) 
fig.savefig('stow_curve_fit.png') 

Answer 2

あなたがフィットしようとしている正弦波の周波数を知っている場合、あなたは正弦波に合わせて線形回帰を使用することができます。任意の正弦波は、正弦関数と余弦関数の線形結合によって表すことができます。線形回帰を使用して正弦と余弦の係数を見つけることができます。このアプローチの良い点は、最初の推測が必要なく、回帰式を満たす回答が1つしかないことです（たとえば、「間違っている」という回答は得られません）。

http://exnumerus.blogspot.com/2010/04/how-to-fit-sine-wave-example-in-python.htmlには、サンプルコードの短いチュートリアルがあります。