Синусоидальная регрессионная линия со сципионом

Question

Я попробовал следующее, чтобы найти синусоидальную регрессию, но не могу нарисовать синусоидальную кривую. Что я тут не так делаю?

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

def sinfunc(x, a, b, c, d):
    return a * np.sin(b * (x - np.radians(c)))+d

year=np.arange(0,24,2)
population=np.array([10.2,11.1,12,11.7,10.6,10,10.6,11.7,12,11.1,10.2,10.2])

popt, pcov = curve_fit(sinfunc, year, population, p0=None)

x_data = np.linspace(0, 25, num=100)

plt.scatter(year,population,label='Population')
plt.plot(x_data, sinfunc(x_data, *popt), 'r-',label='Fitted function')
plt.title("Year vs Population")
plt.xlabel('Year')
plt.ylabel('Population')
plt.legend()
plt.show()

TI-nspire показывает y = sin (0.58x-1) + 11

Обновление

Если я использую p0=[1,0.4,1,5], это работает хорошо. Но не должно ли это быть автоматом c?

Answer 1

То, что вы делаете «неправильно», передает p0=None в curve_fit().

Все методы подбора действительно требуют начальных значений. К сожалению, scipy.optimize.curve_fit() имеет совершенно неоправданную опцию, позволяющую вам не устанавливать начальные значения и молча (даже не предупреждать !!), делая абсурдное предположение, что все значения имеют начальные значения 1,0. Оказывается, что для вашей проблемы эти начальные значения, которые невозможно оправдать и разбить по замыслу, настолько плохи, что не удается найти хороший ответ. Это не редкость. curve_fit лжет вам, что p0=None приемлемо, и вы верите, что l ie.

Решение состоит в том, чтобы признать, что смещение, очевидно, составляет около 11, и использовать p0=[1.0, 0.5, 0.5, 11.0].

Вы можете рассмотреть возможность использования lmfit (https://lmfit.github.io/lmfit-py/). для этой проблемы (отказ от ответственности: я ведущий автор). lmfit имеет класс Model для подгонки кривой, который имеет несколько полезных функций, которые могут быть здесь полезны (не то, что curve_fit не может решить эту проблему - может). С lmfit ваша подгонка может выглядеть следующим образом:

import numpy as np
import matplotlib.pyplot as plt
from lmfit import Model

def sinfunc(x, a, b, c, d):
    return a * np.sin(b*(x - c)) + d

year=np.arange(0,24,2)
population=np.array([10.2,11.1,12,11.7,10.6,10,
                     10.6,11.7,12,11.1,10.2,10.2])

# build model from your model function  
model  = Model(sinfunc)

# create parameters (with initial values!). Note that parameters 
# are named from the argument names of your model function
params = model.make_params(a=1, b=0.5, c=0.5, d=11.0)

# you can set min/max for any parameter to put bounds on the values
params['a'].min = 0
params['c'].min = -np.pi
params['c'].max = np.pi

# do the fit to your data with those parameters
result = model.fit(population, params, x=year)

# print out report of fit statistics and parameter values+uncertainties
print(result.fit_report())

# plot data and fit result
plt.scatter(year,population,label='Population')
plt.plot(year, result.best_fit, 'r-',label='Fitted function')
plt.title("Year vs Population")
plt.xlabel('Year')
plt.ylabel('Population')
plt.legend()
plt.show()

Это распечатает отчет

[[Model]]
    Model(sinfunc)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 26
    # data points      = 12
    # variables        = 4
    chi-square         = 0.00761349
    reduced chi-square = 9.5169e-04
    Akaike info crit   = -80.3528861
    Bayesian info crit = -78.4132595
[[Variables]]
    a:  1.00465520 +/- 0.01247767 (1.24%) (init = 1)
    b:  0.57528444 +/- 0.00198556 (0.35%) (init = 0.5)
    c:  1.80990367 +/- 0.03823815 (2.11%) (init = 0.5)
    d:  11.0250780 +/- 0.00925246 (0.08%) (init = 11)
[[Correlations]] (unreported correlations are < 0.100)
    C(b, c) =  0.812
    C(b, d) =  0.245
    C(c, d) =  0.234

и даст график

Но, опять же, проблема в том, что вы были обмануты, полагая, что p0=None - разумное использование curve_fit().