在Scipy中使用curve_fit擬合矢量函數

Question

我想使用Scipy的curve_fit（或更合適的東西，如果可用）擬合具有矢量輸出的函數。例如，考慮以下功能：

import numpy as np 
def fmodel(x, a, b): 
    return np.vstack([a*np.sin(b*x), a*x**2 - b*x, a*np.exp(b/x)]) 

每個組件都是不同的功能，但它們共享我希望適合的參數。理想情況下，我會做這樣的事情：

x = np.linspace(1, 20, 50) 
a = 0.1 
b = 0.5 
y = fmodel(x, a, b) 
y_noisy = y + 0.2 * np.random.normal(size=y.shape) 

from scipy.optimize import curve_fit 
popt, pcov = curve_fit(f=fmodel, xdata=x, ydata=y_noisy, p0=[0.3, 0.1]) 

但curve_fit不帶矢量輸出功能的工作，和一個錯誤Result from function call is not a proper array of floats.被拋出。我所做的是將這樣的輸出展平：

def fmodel_flat(x, a, b): 
    return fmodel(x[0:len(x)/3], a, b).flatten() 

popt, pcov = curve_fit(f=fmodel_flat, xdata=np.tile(x, 3), 
         ydata=y_noisy.flatten(), p0=[0.3, 0.1]) 

這個工作。如果不是一個矢量函數，我實際上也適用於具有不同輸入的幾個函數，但它們共享模型參數，我可以連接輸入和輸出。

有沒有更合適的方法來適應向量函數與Scipy或可能一些額外的模塊？對我來說主要考慮的是效率 - 實際的功能要複雜得多，而且擬合可能需要一些時間，所以如果這種使用curve_fit被破壞並導致運行時間過長，我想知道我應該做些什麼。

Answer 1

我認爲從效率的角度來看，你所做的是完美無瑕的。我會試着看看實施情況，並提出更多的定量內容，但暫時這裏是我的推理。

你的曲線擬合過程中正在做的是優化參數(a,b)這樣

res = sum_i |f(x_i; a,b)-y_i|^2 

是最小的。我的意思是說你有任意維度的數據點(x_i,y_i)，兩個參數(a,b)和一個近似於查詢點x_i處的數據的擬合模型。

曲線擬合算法從一個起始(a,b)對開始，將其放入一個計算上述平方誤差的黑盒中，並試圖找出一個產生較小誤差的新對(a',b')。我的觀點是，上述誤差對於擬合算法而言實際上是一個黑盒子：配件的配置空間僅由參數(a,b)定義。如果你想象如何實現一個簡單的曲線擬合函數，你可以想象你試圖做一個梯度下降，平方誤差作爲成本函數。

現在，黑箱計算錯誤的擬合程序應該是無關緊要的。很容易看出，x_i的維度對於標量函數來說真的無關緊要，因爲如果您有1000個1d查詢點以適應或者3d空間中的10x10x10網格並不重要。重要的是你有1000點x_i你需要從模型中計算f(x_i) ~ y_i。

應該進一步注意的唯一的微妙之處在於，在矢量值函數的情況下，錯誤的計算並不是微不足道的。在我看來，使用向量值函數的2-範數來定義每個點的誤差是很好的。但是，嘿：在這種情況下，在點x_i的平方誤差是

|f(x_i; a,b)-y_i|^2 == sum_k (f(x_i; a,b)[k]-y_i[k])^2 

這意味着該平方誤差爲每個組件被累積。這意味着你現在正在做的是正確的：通過複製你的點並且單獨考慮函數的每個組成部分，你的平方誤差將精確地包含在每個點的誤差的2-範數。

所以我的觀點是你正在做的是數學上正確的，我不希望擬合過程的任何行爲都依賴於如何處理多元/矢量值函數的方式。

Answer 2

如果我可以如此直言以推薦我自己的包symfit，我認爲它正是你所需要的。在docs中可找到擬合共享參數的示例。上述將成爲

您的具體問題：

from symfit import variables, parameters, Model, Fit, sin, exp 

x, y_1, y_2, y_3 = variables('x, y_1, y_2, y_3') 
a, b = parameters('a, b') 
a.value = 0.3 
b.value = 0.1 

model = Model({ 
    y_1: a * sin(b * x), 
    y_2: a * x**2 - b * x, 
    y_3: a * exp(b/x), 
}) 

xdata = np.linspace(1, 20, 50) 
ydata = model(x=xdata, a=0.1, b=0.5) 
y_noisy = ydata + 0.2 * np.random.normal(size=(len(model), len(xdata))) 

fit = Fit(model, x=xdata, y_1=y_noisy[0], y_2=y_noisy[1], y_3=y_noisy[2]) 
fit_result = fit.execute() 

退房的docs更多！