最小二乘拟合正弦幂级数

Question

我正在尝试拟合以下形式的函数：

其中 A 和 B 是固定常数。在 scipy 中，我通常（并且我认为相当规范）处理此类问题的方法是这样的：

def func(t, coefs):
    phase = np.poly1d(coefs)(t)
    return A * np.cos(phase) + B

def fit(time, data, guess_coefs): 
    residuals = lambda p: func(time, p) - data
    fit_coefs = scipy.optimize.leastsq(residuals, guess_coefs) 
    return fit_coefs

这工作正常，但我想提供一个解析雅可比行列式以改进收敛。因此：

def jacobian(t, coefs):
    phase = np.poly1d(coefs, t)
    the_jacobian = []
    for i in np.arange(len(coefs)):
        the_jac.append(-A*np.sin(phase)*(t**i))
    return the_jac

def fit(time, data, guess_coefs):
    residuals = lambda p: func(time, p) - data
    jac = lambda p: jacobian(time, p)
    fit_coefs = scipy.optimize.leastsq(residuals, guess_coefs, 
                                       Dfun=jac, col_deriv=True)

即使订单数量为 2 或更少，这也不起作用。使用 optimize.check_gradient() 进行快速检查也不会给出肯定的结果。

我几乎可以肯定 Jacobian 和代码是正确的（尽管请纠正我）并且问题更根本：Jacobian 中的 t**i 项导致溢出错误。这不是函数本身的问题，因为这里的单项式项乘以它们的系数，系数很小。

我的问题是：

我上面所做的代码是否有问题？
还是有其他问题？
如果我的假设是正确的，有没有办法对拟合函数进行预处理，使雅可比行列式的表现更好？也许我可以适应数据和时间的对数，或者其他东西。

谢谢！

编辑：忘了原函数形式的正方形

Answer 1

poly1D 函数首先具有最高系数，而您的雅可比函数首先假定最低系数。如果在 jacobian 中你做出 return 声明 return the_jac[::-1]（并修复更明显的拼写错误）你的函数将通过 optimize.check_gradient() 并在 leastsq() 中正常工作。

你关于数值稳定性的进一步问题在这里也是有根据的。如果 t 的值很大且系数很大，则很容易出现数值精度问题：例如，在 32 位精度下，如果多项式的计算结果超过 10^ 8，正弦波的相位将完全丢失在尾数中，正弦评估的结果基本上是垃圾！

您可以通过在多项式中使用模幂来解决这个问题：基本上您关心的多项式的每一项都是 (a_k t ** k) % p，其中 p = 2 * np.pi 是正弦曲线的周期。您可以使用 (a_k * (t % (p / a_k)) ** k) % p 为大型 t 计算更高精度的模指数；对于总体 k 的准确性，事情变得有点复杂。有关此类事情的精彩讨论，请参阅 this answer。

最小二乘拟合正弦幂级数

Least squares fit to sinusoidal power series

python

numpy

scipy

mathematical-optimization

curve-fitting