scipy.optimize.curve_fit无法拟合偏移的偏态高斯曲线。

Question

我正在尝试使用scipy的curve_fit函数拟合一个倾斜和偏移的高斯曲线，但我发现在某些条件下，拟合结果相当糟糕，常常给出接近或正好是一条直线的结果。

下面的代码来自于curve_fit文档。提供的代码是一个任意的测试数据集，但很好地展示了这个问题。

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import math as math
import scipy.special as sp

#def func(x, a, b, c):
#    return a*np.exp(-b*x) + c

def func(x, sigmag, mu, alpha, c,a):
    #normal distribution
    normpdf = (1/(sigmag*np.sqrt(2*math.pi)))*np.exp(-(np.power((x-mu),2)/(2*np.power(sigmag,2))))
    normcdf = (0.5*(1+sp.erf((alpha*((x-mu)/sigmag))/(np.sqrt(2)))))
    return 2*a*normpdf*normcdf + c

x = np.linspace(0,100,100)
y = func(x, 10,30, 0,0,1)
yn = y + 0.001*np.random.normal(size=len(x))

popt, pcov = curve_fit(func, x, yn,) #p0=(9,35,0,9,1))

y_fit= func(x,popt[0],popt[1],popt[2],popt[3],popt[4])

plt.plot(x,yn)
plt.plot(x,y_fit)

问题似乎出现在我将高斯函数从零点（使用mu）移得太远时。我尝试过给出初始值，甚至是与我的原始函数相同的值，但这并不能解决问题。对于mu=10的值，curve_fit可以完美地工作，但如果我使用mu>=30，它就不再适合数据了。

- abradd

Answer 1

你可以使用随机初始值多次调用curve_fit，并选择误差最小的参数。

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import math as math
import scipy.special as sp

def func(x, sigmag, mu, alpha, c,a):
    #normal distribution
    normpdf = (1/(sigmag*np.sqrt(2*math.pi)))*np.exp(-(np.power((x-mu),2)/(2*np.power(sigmag,2))))
    normcdf = (0.5*(1+sp.erf((alpha*((x-mu)/sigmag))/(np.sqrt(2)))))
    return 2*a*normpdf*normcdf + c

x = np.linspace(0,100,100)
y = func(x, 10,30, 0,0,1)
yn = y + 0.001*np.random.normal(size=len(x))

results = []
for i in xrange(50):
    p = np.random.randn(5)*10
    try:
        popt, pcov = curve_fit(func, x, yn, p)
    except:
        pass
    err = np.sum(np.abs(func(x, *popt) - yn))
    results.append((err, popt))
    if err < 0.1:
        break

err, popt = min(results, key=lambda x:x[0])
y_fit= func(x, *popt)

plt.plot(x,yn)
plt.plot(x,y_fit)
print len(results)

- HYRY

Answer 2

提供最小化起点通常会产生奇妙的效果。尝试为最小化器提供有关最大值位置和曲线宽度的信息：

popt, pcov = curve_fit(func, x, yn, p0=(1./np.std(yn), np.argmax(yn) ,0,0,1))

将您代码中的这一行更改为sigma=10和mu=50将产生以下结果

- ev-br

---

好的，现在情况开始变得更加美好了。我觉得我有点高估了curve_fit算法，并对它期望过高。现在我手动估算数值，然后将它们插入到curve_fit中，得到了不错的结果。干杯！ - abradd