Scipy 优化 curve_fit 在拟合自定义函数时为相同的参数给出不同的图

Scipy optimize curve_fit gives different plots for same parameters when fitting custom function

我在 Python 中使用 scipy.optimize 拟合自定义函数时遇到问题,我不知道为什么会这样。我从居中和归一化的二项分布(高斯曲线)生成数据,然后拟合一条曲线。当我在拟合数据上绘制我的函数时,预期的结果就在图片中。但是当我进行拟合时,它失败了。

我确信这是一个 pythonic 的东西,因为它应该给参数 a = 1(这就是我定义它的方式)并且它给了它但是然后拟合不好(见图)。但是,如果我将 sigma 更改为 0.65*sigma in:

p_halfg, p_halfg_cov = optimize.curve_fit(lambda x, a:piecewise_half_gauss(x, a, sigma = 0.65*sigma_fit), x, y, p0=[1])

,它给出了几乎完美的拟合(a 是 5/3,如数学预测的那样)。这些配合应该是相同的,但它们不是! 我在下面给出更多评论。你能告诉我发生了什么事以及问题出在哪里吗?

Plot with a=1 and sigma = sigma_fit

Plot with sigma = 0.65*sigma_fit


我从归一化二项分布生成数据(我可以提供我的代码,但值现在更重要)。它是一个 N = 10 和 p = 0.5 的分布,我将它居中并只取曲线的右侧。然后我用我的半高斯函数拟合它,如果它的参数 a = 1,它应该与二项分布相同(并且西格玛等于分布的西格玛, sqrt(np(1-p)) ) .现在的问题首先是尽管得到了正确的参数a值,但它不适合如图所示的数据。

注意奇怪的东西...如果我设置 sigma = 3* sigma_fit,我得到 a = 1/3 并且非常不合适(低估)。如果我将它设置为 0.2*sigma_fit,我也会得到一个不合适的值,a = 1/0.2 = 5(高估)。等等。为什么? (顺便说一下,a=1/sigma 所以拟合程序应该有效)。

import numpy as np
import matplotlib.pyplot as plt

import math
pi = math.pi

import scipy.optimize as optimize

# define my function
sigma_fit = 1
def piecewise_half_gauss(x, a, sigma=sigma_fit):
    """Half of normal distribution curve, defined as gaussian centered at 0 with constant value of preexponential factor for x < 0
    Arguments: x values as ndarray whose numbers MUST be float type (use linspace or np.arange(start, end, step, dtype=float),
    a as a parameter of width of the distribution,
    sigma being the deviation, second moment
    Returns: Half gaussian curve

    Ex:
    >>> piecewise_half_gauss(5., 1)
    array(0.04839414)

    >>> x = np.linspace(0,10,11)
    ... piecewise_half_gauss(x, 2, 3)
    array([0.06649038, 0.06557329, 0.0628972 , 0.05867755, 0.05324133,
       0.04698531, 0.04032845, 0.03366645, 0.02733501, 0.02158627,
       0.01657952])

    >>> piecewise_half_gauss(np.arange(0,11,1, dtype=float), 1, 2.4)
    array([1.66225950e-01, 1.52405153e-01, 1.17463281e-01, 7.61037856e-02,
       4.14488078e-02, 1.89766470e-02, 7.30345854e-03, 2.36286717e-03,
       6.42616248e-04, 1.46914868e-04, 2.82345875e-05])

    """
    return np.piecewise(x, [x >= 0, x < 0],
                        [lambda x: np.exp(-x ** 2 / (2 * ((a * sigma) ** 2))) / (np.sqrt(2 * pi) * sigma * a),
                         lambda x: 1 / (np.sqrt(2 * pi) * sigma)])


# Create normalized data for binomial distribution Bin(N,p)
n = 10
p = 0.5

x = np.array([0., 1., 2., 3., 4., 5.])
y = np.array([0.25231325, 0.20657662, 0.11337165, 0.0417071 , 0.01028484,
       0.00170007])

# Get the estimate for sigma parameter
sigma_fit = (n*p*(1-p))**0.5

# Get fitting parameters
p_halfg, p_halfg_cov = optimize.curve_fit(lambda x, a:piecewise_half_gauss(x, a, sigma = sigma_fit), x, y, p0=[1])
print(sigma_fit, p_halfg, p_halfg_cov)


## Plot the result
# unpack fitting parameters
a = np.float64(p_halfg)

# unpack uncertainties in fitting parameters from diagonal of covariance matrix
#da = [np.sqrt(p_halfg_cov[j,j]) for j in range(p_halfg.size)] # if we fit more parameters
da = np.float64(np.sqrt(p_halfg_cov[0]))

# create fitting function from fitted parameters
f_fit = np.linspace(0, 10, 50)
y_fit = piecewise_half_gauss(f_fit, a)

# Create figure window to plot data
fig = plt.figure(1, figsize=(10,10))

plt.scatter(x, y, color = 'r', label = 'Original points')
plt.plot(f_fit, y_fit, label = 'Fit')
plt.xlabel('My x values')
plt.ylabel('My y values')

plt.text(5.8, .25, 'a = {0:0.5f}$\pm${1:0.6f}'.format(a, da))
plt.legend()

However, if I plot it manually, it fits EXACTLY!

plt.scatter(x, y, c = 'r', label = 'Original points')
plt.plot(np.linspace(0,5,50), piecewise_half_gauss(np.linspace(0,5,50), 1, sigma_fit), label = 'Fit')
plt.legend()

编辑——已解决: 是绘图问题,需要用到

y_fit = piecewise_half_gauss(f_fit, a, sigma = 0.6*sigma_fit)

问题在于绘制和拟合参数——如果我用不同的 sigma 拟合它,我还需要在生成 y_fit:[=13 时在绘图部分更改它=]

# Get fitting parameters
p_halfg, p_halfg_cov = optimize.curve_fit(lambda x, a:piecewise_half_gauss(x, a, sigma = 0.6*sigma_fit), x, y, p0=[1])
...

y_fit = piecewise_half_gauss(f_fit, a, sigma = 0.6*sigma_fit)