从Scipy-Fmin-Guassian模型到实数d

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. 定义你的损失函数，损失=L（P，x，y），你想在你的参数上最小化
3. 可选：定义一个返回雅可比矩阵的函数，即损失函数的偏导数w.r.t.模型参数。*
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