函数求和的最小二乘法

Least Squares Method for a sum of functions

我想使用 scipy.optimize 模块中的 curve_fit 函数来确定振幅、频率、正弦函数总和(和一个 y0)的相位。当我知道许多要使用的正弦时,这很容易做到。例如,当我从 DFT(离散傅立叶变换)中知道两个频率:1.1520.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

假设您事先知道频率,问题就很简单。您可以将频率的下限设置为 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')