python - Scipy Fmin Guassian 模型到真实数据

Question

我一直在尝试解决这个问题，但实际上还没有看到一个例子或任何我的大脑能够用来前进的东西。

目标是通过最小化实际数据与需要合理估计的未知参数产生的模型之间的总卡方来找到模型高斯曲线(高斯曲线的位置、幅度和宽度未知)。 scipy.optimize.fmin 已经出现，但我以前从未使用过它，而且我对 python 还很陌生......

最终，我想将原始数据与模型一起绘制 - 我之前使用过 pyplot，它只是生成模型并使用 fmin，这让我完全困惑我本质上在这里:

def gaussian(a, b, c, x):
    return a*np.exp(-(x-b)**2/(2*c**2))

我见过多种生成模型的方法，这让我很困惑，因此我没有代码!我已经通过 np.loadtxt 导入了我的数据文件。

感谢任何可以建议框架或提供帮助的人。

Answer 1

像这样的模型拟合问题基本上涉及四个(或五个)主要步骤:

1. 定义您的正向模型，yhat = F(P, x)，它采用一组参数P和自变量x ，并估计您的响应变量y
2. 定义您希望在参数上最小化的损失函数，loss = L(P, x, y)
3. 可选:定义一个返回雅可比矩阵的函数，即损失函数相对于的偏导数。您的模型参数。*
4. 对模型参数进行初步猜测
5. 将所有这些插入优化器之一并获取模型的拟合参数

这是一个可以帮助您入门的示例:

import numpy as np
from scipy.optimize import minimize
from matplotlib import pyplot as pp

# function that defines the model we're fitting
def gaussian(P, x):
    a, b, c = P
    return a*np.exp(-(x-b)**2 /( 2*c**2))

# objective function to minimize
def loss(P, x, y):
    yhat = gaussian(P, x)
    return ((y - yhat)**2).sum()

# generate a gaussian distribution with known parameters
amp = 1.3543
pos = 64.546
var = 12.234
P_real = np.array([amp, pos, var])

# we use the vector of real parameters to generate our fake data
x = np.arange(100)
y = gaussian(P_real, x)
# add some gaussian noise to make things harder
y_noisy = y + np.random.randn(y.size)*0.5

# minimize needs an initial guess at the model parameters
P_guess = np.array([1, 50, 25])

# minimize provides a unified interface to all of scipy's solvers. you
# can also access them individually in scipy.optimize, but the
# standalone versions have annoying differences in their syntax. for now
# we'll use the Nelder-Mead solver, which doesn't use the Jacobian. we
# also need to hand it x and y_noisy as additional args to loss()
res = minimize(loss, P_guess, method='Nelder-Mead', args=(x, y_noisy))

# res is a dict containing the results of the optimization. in particular we
# want the optimized model parameters:
P_fit = res['x']

# we can pass these to gaussian() to evaluate our fitted model
y_fit = gaussian(P_fit, x)

# now let's plot the results:
fig, ax = pp.subplots(1,1)
ax.hold(True)
ax.plot(x, y, '-r', lw=2, label='Real')
ax.plot(x, y_noisy, '-k', alpha=0.5, label='Noisy')
ax.plot(x, y_fit, '--b', lw=5, label='Fit')
ax.legend(loc=0, fancybox=True)

*一些求解器，例如共轭梯度方法，将雅可比行列式作为附加参数，总的来说，这些求解器更快、更稳健，但如果您感到懒惰并且性能不是那么重要，那么您通常可以在不提供雅可比行列式的情况下逃脱，在这种情况下，它将使用有限差分法来估计梯度。

您可以阅读有关不同求解器的更多信息 here