python - 函数和的最小二乘法

Question

我想使用 curve_fit函数来自 scipy.optimize模块来确定正弦函数之和的幅度、频率、相位(以及一个 y0 )。当我知道要使用的正弦数时，这很容易做到。例如，当我从 DFT(离散傅里叶变换)知道两个频率时:1.152和0.432我可以定义一个函数:

def func(x, amp1, amp2, freq1 , freq2, phase1, phase2, y0):
    return amp1*np.sin(freq1*x + phase1) + amp2*np.sin(freq2*x + phase2) + y0

然后，使用curve_fit并限制频率间隔我可以找到一个很好的拟合:

param, _ = curve_fit(func, t, data, bounds=([-np.inf, -np.inf, 1.14, 0.43, -np.inf, -np.inf, -np.inf], [np.inf, np.inf, 1.16, 0.44, np.inf, np.inf, np.inf]))

看起来很棒:

但在这种情况下，我已经准备好了数据并且我已经知道了一些频率。你知道如何定义func吗？仅一次并处理所有情况(例如五个正弦函数)？我尝试将参数放入列表中，例如amp = [amp1, amp2, ... ]我已经迭代了它们的长度。但是定义bounds有一个问题对于参数列表。 bounds确保现实模型非常重要。

解决方案不必基于curve_fit .

Answer 1

假设您事先知道频率，问题就很简单。您可以将频率的下限设置为 0，将上限设置为 2 * pi * freq。对于安培，设置任意数字(如果不需要边界，则设置 np.inf)。

您可以将函数表述为lambda x, amp1,phase1, amp2,phase2...: y，curve_fit可以接受未定义数量的函数只要您提供正确的初始猜测即可进行论证。

五个频率的示例代码:

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

x = np.linspace(0,10,60)
w = [1,2,3,4,5]

a = [1,4,2,3,0.1]
x0 = [0,1,0,1,0.5]

y = np.sum(a_i * np.sin(w_i * x - x0_i) for w_i, a_i, x0_i in zip(w,a, x0))  #base_data
yr = y + np.random.normal(0,0.5, size=x.size)   #noisy data

def func(x, *args):
    """ function of the form lambda x, amp1, phase1, amp2, phase2...."""
    return np.sum(a_i * np.sin(w_i * (x-x0)) for w_i, a_i, x0
                  in zip(w,args[::2], args[1::2]))

ubounds = np.zeros(len(w) * 2)
ubounds[::2] = 10   #setting amp max value to 10 (arbitrary)
ubounds[1::2] = np.asarray(w) * 2 * np.pi
p0 = [0] * 10   # note p0 size
popt, pcov = curve_fit(func, x, yr, p0, bounds=(0, ubounds))
amps, phases = popt[::2], popt[1::2]

plt.plot(x,func(x, *popt))
plt.plot(x,yr, 'go')