使用 LeastSq 将 CDF 与 2 高斯拟合答案

【问题标题】：Fit CDF with 2 Gaussian using LeastSq使用 LeastSq 将 CDF 与 2 高斯拟合
【发布时间】：2018-08-28 16:09:58
【问题描述】：

我正在尝试将经验 CDF 图拟合到两个高斯 cdf，因为它似乎有两个峰值，但它不起作用。我使用 scipy.optimize 中的 leastsq 和 scipy.special 中的 erf 函数拟合曲线。拟合仅给出值为 2 的常数线。我不确定我在代码的哪个部分犯了错误。任何指针都会有所帮助。谢谢！

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

x = np.array([ 90.64115156,  90.85690063,  91.07264971,  91.28839878,
        91.50414786,  91.71989693,  91.93564601,  92.15139508,
        92.36714415,  92.58289323,  92.7986423 ,  93.01439138,
        93.23014045,  93.44588953,  93.6616386 ,  93.87738768,
        94.09313675,  94.30888582,  94.5246349 ,  94.74038397,
        94.95613305,  95.17188212,  95.3876312 ,  95.60338027,
        95.81912935,  96.03487842,  96.2506275 ,  96.46637657,
        96.68212564,  96.89787472,  97.11362379,  97.32937287,
        97.54512194,  97.76087102,  97.97662009,  98.19236917,
        98.40811824,  98.62386731,  98.83961639,  99.05536546,
        99.27111454,  99.48686361,  99.70261269,  99.91836176,
       100.13411084, 100.34985991, 100.56560899, 100.78135806,
       100.99710713, 101.21285621])
y = np.array([3.33333333e-04, 3.33333333e-04, 3.33333333e-04, 1.00000000e-03,
       1.33333333e-03, 3.33333333e-03, 6.66666667e-03, 1.30000000e-02,
       2.36666667e-02, 3.40000000e-02, 5.13333333e-02, 7.36666667e-02,
       1.01666667e-01, 1.38666667e-01, 2.14000000e-01, 3.31000000e-01,
       4.49666667e-01, 5.50000000e-01, 6.09000000e-01, 6.36000000e-01,
       6.47000000e-01, 6.54666667e-01, 6.61000000e-01, 6.67000000e-01,
       6.76333333e-01, 6.84000000e-01, 6.95666667e-01, 7.10000000e-01,
       7.27666667e-01, 7.50666667e-01, 7.75333333e-01, 7.93333333e-01,
       8.11333333e-01, 8.31333333e-01, 8.56333333e-01, 8.81333333e-01,
       9.00666667e-01, 9.22666667e-01, 9.37666667e-01, 9.47333333e-01,
       9.59000000e-01, 9.70333333e-01, 9.77333333e-01, 9.83333333e-01,
       9.90333333e-01, 9.93666667e-01, 9.96333333e-01, 9.99000000e-01,
       9.99666667e-01, 1.00000000e+00])
plt.plot(a,b,'r.')

# Fitting with 2 Gaussian
from scipy.special import erf
from scipy.optimize import leastsq

def two_gaussian_cdf(params, x):
    (mu1, sigma1, mu2, sigma2) = params
    model = 0.5*(1 + erf( (x-mu1)/(sigma1*np.sqrt(2)) )) +\
            0.5*(1 + erf( (x-mu2)/(sigma2*np.sqrt(2)) ))
    return model

def residual_two_gaussian_cdf(params, x, y):
    model = double_gaussian(params, x)
    return model - y

params = [5.,2.,1.,2.]
out = leastsq(residual_two_gaussian_cdf,params,args=(x,y))
double_gaussian(out[0],x)
plt.plot(x,two_gaussian_cdf(out[0],x))

回到这个情节

【问题讨论】：

至少部分问题是初始参数值。尝试“params = [100.,2.,100.,2.]”，它会更现实一些。我不知道是否还有其他问题，但是您帖子中的原始值肯定会给您带来麻烦。我使用 scipy 模块 scipy.optimize.differential_evolution 来搜索初始参数，这些是我从该模块收到的结果的整数。
哦，是的。我只是意识到 x 轴值是如此之大，以至于 1 数量级的参数不能很好地工作。我尝试了您提供的值并且它有效。很高兴知道 scipy.optimize.differential_evolution 模块。谢谢！
@JamesPhillips，你能告诉我如何运行 scipy.optimize.differential_evolution 模块来获得估计吗？
在下面查看我的答案。我无法在 cmets 中格式化代码，因此以答案的形式提供了代码。

标签： python-3.x curve-fitting gaussian least-squares cdf

【解决方案1】：

您可能会发现lmfit（请参阅http://lmfit.github.io/lmfit-py/）是此处leastsq 的有用替代品，因为它为优化和曲线拟合提供了更高级别的接口（尽管仍然基于scipy.optimize.leastsq）。使用 lmfit，您的示例可能如下所示（删除 x 和 y 数据的定义）：

#!/usr/bin/env python  
import numpy as np
from scipy.special import erf
import matplotlib.pyplot as plt
from lmfit import Model

# define the basic model.  I included an amplitude parameter
def gaussian_cdf(x, amp, mu, sigma):
    return (amp/2.0)*(1 + erf( (x-mu)/(sigma*np.sqrt(2))))

# create a model that is the sum of two gaussian_cdfs
# note that a prefix names each component and will be
# applied to the parameter names for each model component
model = Model(gaussian_cdf, prefix='g1_') + Model(gaussian_cdf, prefix='g2_')

# make a parameters object -- a dict with parameter names
# taken from the arguments of your model function and prefix
params = model.make_params(g1_amp=0.50, g1_mu=94, g1_sigma=1,
                           g2_amp=0.50, g2_mu=98, g2_sigma=1.)

# you can apply bounds to any parameter
#params['g1_sigma'].min = 0   # sigma must be > 0!

# you may want to fix the amplitudes to 0.5:
#params['g1_amp'].vary = False
#params['g2_amp'].vary = False

# run the fit
result = model.fit(y, params, x=x)

# print results
print(result.fit_report())

# plot results, including individual components
comps = result.eval_components(result.params, x=x)

plt.plot(x, y,'r.', label='data')
plt.plot(x, result.best_fit, 'k-', label='fit')
plt.plot(x, comps['g1_'], 'b--', label='g1_')
plt.plot(x, comps['g2_'], 'g--', label='g2_')
plt.legend()
plt.show()

这会打印出一份报告

[[Model]]
    (Model(gaussian_cdf, prefix='g1_') + Model(gaussian_cdf, prefix='g2_'))
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 66
    # data points      = 50
    # variables        = 6
    chi-square         = 0.00626332
    reduced chi-square = 1.4235e-04
    Akaike info crit   = -437.253376
    Bayesian info crit = -425.781238
[[Variables]]
    g1_amp:    0.65818908 +/- 0.00851338 (1.29%) (init = 0.5)
    g1_mu:     93.8438526 +/- 0.01623273 (0.02%) (init = 94)
    g1_sigma:  0.54362156 +/- 0.02021614 (3.72%) (init = 1)
    g2_amp:    0.34058664 +/- 0.01153346 (3.39%) (init = 0.5)
    g2_mu:     97.7056728 +/- 0.06408910 (0.07%) (init = 98)
    g2_sigma:  1.24891832 +/- 0.09204020 (7.37%) (init = 1)
[[Correlations]] (unreported correlations are < 0.100)
    C(g1_amp, g2_amp)     = -0.892
    C(g2_amp, g2_sigma)   =  0.848
    C(g1_amp, g2_sigma)   = -0.744
    C(g1_amp, g1_mu)      =  0.692
    C(g1_amp, g2_mu)      =  0.662
    C(g1_mu, g2_amp)      = -0.607
    C(g1_amp, g1_sigma)   =  0.571

还有这样的情节：

这种搭配并不完美，但它应该可以帮助您入门。

【讨论】：

这是一个很好的答案！我一定会研究lmfit。该报告提供了配件的详细信息。谢谢！我很感激。

【解决方案2】：

以下是我如何使用 scipy.optimize.differential_evolution 模块为曲线拟合生成初始参数估计值。我已将误差平方和编码为遗传算法的目标，如下所示。这个 scipy 模块使用拉丁超立方算法来确保对参数空间的彻底搜索，这需要在其中搜索的参数范围。在这种情况下，参数边界是自动从数据中导出的，因此无需在代码中手动提供。

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import warnings

from scipy.optimize import differential_evolution
from scipy.special import erf


# bounds on parameters are set in generate_Initial_Parameters() below
def two_gaussian_cdf(x, mu1, sigma1, mu2, sigma2):
    model = 0.5*(1 + erf( (x-mu1)/(sigma1*np.sqrt(2)) )) +\
            0.5*(1 + erf( (x-mu2)/(sigma2*np.sqrt(2)) ))
    return model


# function for genetic algorithm to minimize (sum of squared error)
# bounds on parameters are set in generate_Initial_Parameters() below
def sumOfSquaredError(parameterTuple):
    warnings.filterwarnings("ignore") # do not print warnings by genetic algorithm
    return np.sum((yData - two_gaussian_cdf(xData, *parameterTuple)) ** 2)


def generate_Initial_Parameters():
    # data min and max used for bounds
    maxX = max(xData)
    minX = min(xData)
    maxY = max(yData)
    minY = min(yData)

    parameterBounds = []
    parameterBounds.append([minX, maxX]) # parameter bounds for mu1
    parameterBounds.append([minY, maxY]) # parameter bounds for sigma1
    parameterBounds.append([minX, maxX]) # parameter bounds for mu2
    parameterBounds.append([minY, maxY]) # parameter bounds for sigma2

    # "seed" the numpy random number generator for repeatable results
    result = differential_evolution(sumOfSquaredError, parameterBounds, seed=3)
    return result.x


xData = np.array([ 90.64115156,  90.85690063,  91.07264971,  91.28839878,
        91.50414786,  91.71989693,  91.93564601,  92.15139508,
        92.36714415,  92.58289323,  92.7986423 ,  93.01439138,
        93.23014045,  93.44588953,  93.6616386 ,  93.87738768,
        94.09313675,  94.30888582,  94.5246349 ,  94.74038397,
        94.95613305,  95.17188212,  95.3876312 ,  95.60338027,
        95.81912935,  96.03487842,  96.2506275 ,  96.46637657,
        96.68212564,  96.89787472,  97.11362379,  97.32937287,
        97.54512194,  97.76087102,  97.97662009,  98.19236917,
        98.40811824,  98.62386731,  98.83961639,  99.05536546,
        99.27111454,  99.48686361,  99.70261269,  99.91836176,
       100.13411084, 100.34985991, 100.56560899, 100.78135806,
       100.99710713, 101.21285621])
yData = np.array([3.33333333e-04, 3.33333333e-04, 3.33333333e-04, 1.00000000e-03,
       1.33333333e-03, 3.33333333e-03, 6.66666667e-03, 1.30000000e-02,
       2.36666667e-02, 3.40000000e-02, 5.13333333e-02, 7.36666667e-02,
       1.01666667e-01, 1.38666667e-01, 2.14000000e-01, 3.31000000e-01,
       4.49666667e-01, 5.50000000e-01, 6.09000000e-01, 6.36000000e-01,
       6.47000000e-01, 6.54666667e-01, 6.61000000e-01, 6.67000000e-01,
       6.76333333e-01, 6.84000000e-01, 6.95666667e-01, 7.10000000e-01,
       7.27666667e-01, 7.50666667e-01, 7.75333333e-01, 7.93333333e-01,
       8.11333333e-01, 8.31333333e-01, 8.56333333e-01, 8.81333333e-01,
       9.00666667e-01, 9.22666667e-01, 9.37666667e-01, 9.47333333e-01,
       9.59000000e-01, 9.70333333e-01, 9.77333333e-01, 9.83333333e-01,
       9.90333333e-01, 9.93666667e-01, 9.96333333e-01, 9.99000000e-01,
       9.99666667e-01, 1.00000000e+00])

# generate initial parameter values
initialParameters = generate_Initial_Parameters()

# curve fit the data
fittedParameters, niepewnosci = curve_fit(two_gaussian_cdf, xData, yData, initialParameters)

# create values for display of fitted peak function
mu1, sigma1, mu2, sigma2 = fittedParameters
y_fit = two_gaussian_cdf(xData, mu1, sigma1, mu2, sigma2)

plt.plot(xData, yData) # plot the raw data
plt.plot(xData, y_fit) # plot the equation using the fitted parameters
plt.show()

print(fittedParameters)

【讨论】：

感谢您撰写解释。我现在明白我们需要将目标函数设置为模型和实际数据之间的差异。谢谢！
如果您有一个大型数据集，您可以在数据的代表性子集上生成初始参数估计，并且与使用整个数据集相比具有更好的性能。我在我的开源拟合库中使用了这种技术。
我只是看看你的图书馆。这是一个庞大的曲线拟合库。很有帮助。
谢谢。您可能有兴趣查看由该库提供支持的在线曲线拟合和曲面拟合 Web 应用程序，URL 是 zunzun.com - 它也是开源的。该库本身带有许多示例，包括用于提高性能的并行处理。