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

あなたがしていることは、効率の立場からはまったく問題ないと思います。私は実装を見てより定量的なものを考案しようとしますが、当分の間は私の推論です。

あなたは

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

が最小であることを(a,b)なパラメータを最適化されたカーブフィッティングの際に何をやっています。これは、任意の次元のデータポイント(x_i,y_i)と、2つのパラメータ(a,b)と、クエリポイントx_iのデータを近似するフィッティングモデルがあることを意味します。

カーブフィッティングアルゴリズムは、最初の(a,b)ペアから開始し、これを上記の2乗誤差を計算するブラックボックスに入れて、小さな誤差を生じる新しい(a',b')ペアを生成しようとします。私の指摘は、上記のエラーは実際にフィッティングアルゴリズムのブラックボックスであるということです。フィッティングのコンフィグレーション空間は、単に(a,b)パラメータで定義されています。あなたが単純なカーブフィッティング関数をどのように実装するか想像すれば、コスト関数としての二乗誤差を使って、例えば勾配降下を試みようと想像することができます。

ここで、ブラックボックスがエラーをどのように計算するかについてのフィッティング手順とは関係ありません。 x_iの次元数はスカラー関数にとっては無関係であることは容易に分かります。なぜなら、1000個の1次元クエリポイントがあるかどうかは関係ありませんし、3次元空間内の10x10x10グリッドも問題ではないからです。重要なのは、モデルからf(x_i) ~ y_iを計算する必要がある1000点のx_iがあることです。

さらに注意する必要があるのは、ベクトル値関数の場合、誤差の計算は簡単ではないということです。私の意見では、ベクトル値関数の2ノルムを使って、それぞれx_i点に誤差を定義するのは良いことです。ちょっとこの場合、点x_iにおける二乗誤差は、各成分の二乗誤差が蓄積されることを意味

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

あります。これはちょうどあなたが今行っていることがちょうど正しいことを意味します：あなたのx_iポイントを複製し、関数の各コンポーネントを個別に考慮すると、二乗誤差は正確に各点の誤差の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をチェックしてください！