Làm cách nào để tôi vừa vặn với đường cong sin vào dữ liệu của tôi bằng giá treo tường và có khối u?

Question

Đối với một dự án trường học, tôi đang cố gắng chứng minh rằng các nền kinh tế theo một mô hình tăng trưởng tương đối hình sin. Ngoài kinh tế học của nó, được thừa nhận là tinh ranh, tôi đang xây dựng một mô phỏng python để chỉ ra rằng ngay cả khi chúng ta để một mức độ ngẫu nhiên nào đó nắm giữ, chúng ta vẫn có thể tạo ra thứ gì đó tương đối hình sin. Tôi hài lòng với dữ liệu của tôi mà tôi đang sản xuất nhưng bây giờ, tôi muốn tìm một cách nào đó để có được một biểu đồ sin phù hợp chặt chẽ với dữ liệu. Tôi biết bạn có thể làm phù hợp đa thức, nhưng bạn có thể làm phù hợp với sin?

Cảm ơn sự giúp đỡ của bạn trước. Hãy cho tôi biết nếu có bất kỳ phần nào của mã bạn muốn xem.

Answer 1

Bạn có thể sử dụng chức năng least-square optimization trong scipy để phù hợp với bất kỳ chức năng tùy ý nào khác. Trong trường hợp lắp đặt một hàm sin, 3 thông số phù hợp là offset ('a'), biên độ ('b') và pha ('c'). Miễn là bạn cung cấp một dự đoán hợp lý của các thông số, tối ưu hóa nên hội tụ tốt.May mắn thay cho một hàm sin, ước tính đầu tiên của 2 trong số này là dễ dàng: bù đắp có thể được ước tính bằng cách lấy trung bình của dữ liệu và biên độ thông qua RMS (3 * độ lệch chuẩn/sqrt (2)).

Điều này dẫn đến đoạn mã sau:

import numpy as np 
from scipy.optimize import leastsq 
import pylab as plt 

N = 1000 # number of data points 
t = np.linspace(0, 4*np.pi, N) 
data = 3.0*np.sin(t+0.001) + 0.5 + np.random.randn(N) # create artificial data with noise 

guess_mean = np.mean(data) 
guess_std = 3*np.std(data)/(2**0.5) 
guess_phase = 0 

# we'll use this to plot our first estimate. This might already be good enough for you 
data_first_guess = guess_std*np.sin(t+guess_phase) + guess_mean 

# Define the function to optimize, in this case, we want to minimize the difference 
# between the actual data and our "guessed" parameters 
optimize_func = lambda x: x[0]*np.sin(t+x[1]) + x[2] - data 
est_std, est_phase, est_mean = leastsq(optimize_func, [guess_std, guess_phase, guess_mean])[0] 

# recreate the fitted curve using the optimized parameters 
data_fit = est_std*np.sin(t+est_phase) + est_mean 

plt.plot(data, '.') 
plt.plot(data_fit, label='after fitting') 
plt.plot(data_first_guess, label='first guess') 
plt.legend() 
plt.show() 

Edit: Tôi cho rằng bạn biết số kỳ trong sóng sin. Nếu bạn không, nó hơi phức tạp hơn để phù hợp. Bạn có thể thử và đoán số thời gian bằng cách vẽ theo cách thủ công và thử và tối ưu hóa nó theo thông số thứ 4 của bạn.

Answer 2

Các phương pháp hiện tại để phù hợp với đường cong tội lỗi đối với tập dữ liệu đã cho yêu cầu dự đoán đầu tiên về các tham số, theo sau là một quá trình tương tác. Đây là một vấn đề hồi quy phi tuyến tính.

Một phương pháp khác bao gồm việc chuyển đổi hồi quy phi tuyến tính thành hồi quy tuyến tính nhờ phương trình tích phân thuận tiện. Sau đó, không cần đoán ban đầu và không cần quá trình lặp: khớp nối được lấy trực tiếp.

Trong trường hợp của hàm y = a + r*sin(w*x+phi) hoặc y=a+b*sin(w*x)+c*cos(w*x), xem các trang 35-36 của giấy "Régression sinusoidale" công bố trên Scribd

Trong trường hợp của hàm y = a + p*x + r*sin(w*x+phi): trang 49-51 của chương "tuyến tính hỗn hợp và sin hồi quy" .

Trong trường hợp các chức năng phức tạp hơn, quá trình chung được giải thích trong chương "Generalized sinusoidal regression" trang 54-61, theo sau là một ví dụ bằng số y = r*sin(w*x+phi)+(b/x)+c*ln(x), các trang 62-63

Answer 3

More userfriendly đối với chúng tôi là chức năng curvefit. Ở đây một ví dụ:

import numpy as np 
from scipy.optimize import curve_fit 
import pylab as plt 

N = 1000 # number of data points 
t = np.linspace(0, 4*np.pi, N) 
data = 3.0*np.sin(t+0.001) + 0.5 + np.random.randn(N) # create artificial data with noise 

guess_freq = 1 
guess_amplitude = 3*np.std(data)/(2**0.5) 
guess_phase = 0 
guess_offset = np.mean(data) 

p0=[guess_freq, guess_amplitude, 
    guess_phase, guess_offset] 

# create the function we want to fit 
def my_sin(x, freq, amplitude, phase, offset): 
    return np.sin(x * freq + phase) * amplitude + offset 

# now do the fit 
fit = curve_fit(my_sin, t, data, p0=p0) 

# we'll use this to plot our first estimate. This might already be good enough for you 
data_first_guess = my_sin(t, *p0) 

# recreate the fitted curve using the optimized parameters 
data_fit = my_sin(t, *fit[0]) 

plt.plot(data, '.') 
plt.plot(data_fit, label='after fitting') 
plt.plot(data_first_guess, label='first guess') 
plt.legend() 
plt.show() 

Answer 4

Đây là một chức năng phù hợp tham số miễn fit_sin() mà không yêu cầu thủ đoán của tần số:

import numpy, scipy.optimize 

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)} 

Các dự đoán tần số ban đầu được đưa ra bởi tần số đỉnh trong miền tần số sử dụng FFT. Kết quả phù hợp là giả thiết gần như hoàn hảo chỉ có một tần số thống trị (trừ tần số không tần số).

import pylab as plt 

N, amp, omega, phase, offset, noise = 500, 1., 2., .5, 4., 3 
#N, amp, omega, phase, offset, noise = 50, 1., .4, .5, 4., .2 
#N, amp, omega, phase, offset, noise = 200, 1., 20, .5, 4., 1 
tt = numpy.linspace(0, 10, N) 
tt2 = numpy.linspace(0, 10, 10*N) 
yy = amp*numpy.sin(omega*tt + phase) + offset 
yynoise = yy + noise*(numpy.random.random(len(tt))-0.5) 

res = fit_sin(tt, yynoise) 
print("Amplitude=%(amp)s, Angular freq.=%(omega)s, phase=%(phase)s, offset=%(offset)s, Max. Cov.=%(maxcov)s" % res) 

plt.plot(tt, yy, "-k", label="y", linewidth=2) 
plt.plot(tt, yynoise, "ok", label="y with noise") 
plt.plot(tt2, res["fitfunc"](tt2), "r-", label="y fit curve", linewidth=2) 
plt.legend(loc="best") 
plt.show() 

Kết quả là tốt ngay cả với tiếng ồn cao:

Amplitude = 1,00660540618, góc freq = 2,03370472482, giai đoạn = ,360276844224, bù đắp = 3,95747467506, Max.. Cov. = 0,0122923578658