適合正弦電源系列

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（）的快速檢查也沒有給出肯定的結果。

我幾乎可以確定雅可比矩陣和代碼是正確的（儘管請糾正我），而且問題更爲根本：雅可比矩陣中的t ** i項導致溢出錯誤。這在函數本身中不是問題，因爲這裏單項項乘以它們的非常小的係數。

我的問題是：

1. 有什麼錯誤的代碼明智與我之前所做的那樣？
2. 或者還有其他一些問題嗎？
3. 如果我的假設是正確的，是否有辦法預適應函數，這樣雅可比行爲更好？也許我可以適應數據和時間的對數，或者什麼的。

謝謝！

編輯：忘了原來的功能形式

Answer 1

的poly1D功能具有最高係數第一方，而你的雅可比函數首先假定係數最低。如果在jacobian中，您將返回語句return the_jac[::-1]（並修復了更明顯的拼寫錯誤），您的功能將通過optimize.check_gradient()並在leastsq()中正常工作。

關於數值穩定性的進一步問題也在這裏保證。如果你有大的值和大量的係數，你可能很容易出現數值精度問題：例如，如果你的多項式的計算結果超過10^8，那麼在32位精度下，正弦曲線的相位將會完全丟失在尾數和正弦評估的結果將基本上是垃圾！

你可以通過使用多項式內模冪解決這個問題：基本上所有你關心的多項式的每屆任期爲(a_k t ** k) % p，其中p = 2 * np.pi是正弦曲線的週期。您可以使用(a_k * (t % (p/a_k)) ** k) % p來計算此模塊指數，從而獲得較高的精度，適用於大型t;對於精確度大的k，情況會變得更復雜一些。請參閱this answer以獲得有關此類事情的良好討論。