最小二乘拟合正弦幂级数

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和代码是正确的(虽然请纠正我)并且问题更为基础:雅各比派中的t**i术语会导致溢出错误.这在函数本身中不是问题,因为这里单项式项乘以它们非常小的系数.

我的问题是:

1. 我上面做了什么,代码有什么不对吗？
2. 还是有其他问题吗？
3. 如果我的假设是正确的,有没有办法预处理拟合函数,那么雅可比行为更好？也许我可以适应数据和时间的对数,或者其他东西.

谢谢!

编辑:忘记原始功能形式的方块

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正弦波的周期在哪里.你可以计算这个模指数为大型高精密t带(a_k * (t % (p / a_k)) ** k) % p; 对于整体的准确性k,事情变得有点复杂.看到这个答案,对这些事情进行了很好的讨论.