Scipyのcurve_fitを用いた正弦波のフィッティングがうまくできない

Question

実現したいこと

y = a sin(bx + c) + d という関数でフィッティングをしたいのですが、scipyのcurve_fitでうまくできませんでした。改善点、アドバイスいただけると幸いです。

発生している問題・エラーメッセージ

現在の出力
[-0.24459067 0.88349125 1.9100531 5.00003073]

正しい出力
[2 3 4 5 ]

該当のソースコード

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# サンプルデータ
x = np.arange(0,10, 0.1)
y = 2 * np.sin(3 * x + 4) + 5

# フィッティング関数
def func(x, a, b, c, d):
    return a * np.sin(b * x + c) + d

# フィッティングを実行
param, cov = curve_fit(func, x, y)
print(param)
y_fit = func(x, param[0], param[1], param[2], param[3])

# グラフの描画
plt.scatter(x, y, c="r", s=5, label="data")
plt.plot(x, y_fit, c='b', linewidth=1, label="fitting")
plt.legend();

Answer 1

天下り的にはなりますが，初期値p0を角速度bに近い値で始めさせればうまくいきます．
適当にa,b,c,d全て3.5からスタートするようにするとフィットしました．

param, cov = curve_fit(func, x, y, p0 = [3.5] * 4)

初期値を乱数にして共分散の小さいものを選ぶのも良いかもしれませんね．

import numpy as np
min_cov = 1e8
retry, max_retry = 0, 100
while min_cov > 1e-10 and retry < max_retry:
    print("try init prams:", p0 := np.random.randint(1, 10, size = 4))
    param, cov = curve_fit(func, x, y, p0 = p0)
    min_cov = min(np.abs(cov).sum(), min_cov)
    retry += 1

if retry == max_retry:
    print("max retried error")
    exit()

print("retried", retry, "times")
print("fitted params:", param)
y_fit = func(x, *param)

最小二乗法でうまくいかない理由

元の式をy = 2 * np.sin(-3 * x + 4) + 5にして乱数[1, 10]の範囲による初期値選択を行うとその片鱗が見えてきます．このときの解は正しくは[2, -3, 4, 5]ですが，-3は乱数の範囲にないので代わりに3が選ばれ，他のパラメータに皺寄せがいく形になり[-2, 3, 2.28, 5]となりました．

今回うまくいかなかった理由は，三角関数自体の周期性などによって，「最小二乗法で得られる誤差関数が描く超平面の局所的最適解が多すぎる．」というのが考えられます．質問者はその局所的最適解に陥ったものと思われます．また，大域的最適解も複数存在しており，先述のような別解を得ることになっているものと思われます．

Comments

[cream_puff]

ご回答ありがとうございました。「最小二乗法でうまくいかない理由」納得がいきました。

Answer 2

curve_fit failing on even a sine waveにも記載されているとおり、最小二乗法をもとにしたcurve_fitは（私は理論を含めた理由は理解できていませんが）正弦関数のフィッティングには向いていません。
よって代わりにFFTをベースにした方法でフィッティングするとよいでしょう。
以下はHow do I fit a sine curve to my data with pylab and numpy?の回答をもとにしたコード例です。

import numpy as np
import matplotlib.pyplot as plt
#from scipy.optimize import curve_fit
import numpy, scipy.optimize


# サンプルデータ
x = np.arange(0,10, 0.1)
y = 2 * np.sin(3 * x + 4) + 5


def fit_sin(tt, yy):
    '''Fit sin to the input time sequence, and return fitting parameters "amp", "omega", "phase", "offset", "freq", "period" and "fitfunc"'''
    tt = numpy.array(tt)
    yy = numpy.array(yy)
    ff = numpy.fft.fftfreq(len(tt), (tt[1]-tt[0]))   # assume uniform spacing
    Fyy = abs(numpy.fft.fft(yy))
    guess_freq = abs(ff[numpy.argmax(Fyy[1:])+1])   # excluding the zero frequency "peak", which is related to offset
    guess_amp = numpy.std(yy) * 2.**0.5
    guess_offset = numpy.mean(yy)
    guess = numpy.array([guess_amp, 2.*numpy.pi*guess_freq, 0., guess_offset])

    def sinfunc(t, A, w, p, c):  return A * numpy.sin(w*t + p) + c
    popt, pcov = scipy.optimize.curve_fit(sinfunc, tt, yy, p0=guess)
    A, w, p, c = popt
    f = w/(2.*numpy.pi)
    fitfunc = lambda t: A * numpy.sin(w*t + p) + c
    return {"amp": A, "omega": w, "phase": p, "offset": c, "freq": f, "period": 1./f, "fitfunc": fitfunc, "maxcov": numpy.max(pcov), "rawres": (guess,popt,pcov)}


res = fit_sin(x, y)
y_fit = res['fitfunc'](x)


# グラフの描画
plt.scatter(x, y, c="r", s=5, label="data")
plt.plot(x, y_fit, c='b', linewidth=1, label="fitting")
plt.legend()
plt.show()

Comments

[cream_puff]

ご回答ありがとうございました。初期値の与え方が非常に参考になりました。

Answer 3

scipyのcurve_fitでは、パラメータをある程度推測して初期値として設定しておくことが重要です。
デフォルトの初期値はすべて1になっているので、そこから正解があまりにも離れすぎていると、うまくフィッティングできません。
使う関数の性質を考慮したり、そのパラメータを使って実際に描画したりして、大体合っている場合にうまくいきやすいです。

以下、初期値を正解±0.3で設定した例です。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# サンプルデータ
x = np.arange(0, 10, 0.1)
y = 2 * np.sin(3 * x + 4) + 5


# フィッティング関数
def func(x, a, b, c, d):
    return a * np.sin(b * x + c) + d


# 初期値
p0 = [2.3, 2.7, 4.3, 4.7]


# フィッティングを実行
param, cov = curve_fit(func, x, y, p0=p0)  # 初期値を設定
print(param)
# y_fit = func(x, param[0], param[1], param[2], param[3])
y_fit = func(x, *param)  # tips:アンパックが使える

# グラフの描画
plt.scatter(x, y, c="r", s=5, label="data")
plt.plot(x, y_fit, c="b", linewidth=1, label="fitting")
plt.legend()

出力: [2. 3. 4. 5.]

グラフ:

（参考）初期値のプロットも重ねた場合:

これくらいズレてしまっていても、うまくフィッティングできる場合もあるようです。

Answer 4

どんな場合にでも使えるものではないですが、ある程度範囲を絞ってフィッティングして、それを初期値に全体をフィッティングするという方法もあります。

p0 = [1, 1, 1, 1]
param0, _ = curve_fit(func, x[:30], y[:30], p0=p0)
param, cov = curve_fit(func, x, y, p0=param0)

(今回の例では param0 の時点で全体にもフィッティングできています。)

このやり方でも、ある程度は初期値(p0)が正解に近くないと、正しく収束できません。
たとえば、p0 = [1, 8, 1, 1] ではうまくいきません。