Установка векторной функции с curve_fit в Scipy

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 для соответствия или сетка 10x10x10 в 3d пространстве. Важно то, что у вас есть 1000 очков x_i, для которых вам необходимо вычислить f(x_i) ~ y_i.

Единственная тонкость, которую следует дополнительно отметить, заключается в том, что в случае вектор-функции вычисление ошибки не является тривиальным. На мой взгляд, хорошо определить ошибку на каждой точке x_i, используя 2-норму вектор-функции.Но погодите: в этом случае погрешность в точке 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 больше!