自定义函数拟合时,Scipy optimize curve_fit对相同参数给出了不同的曲线图。

问题描述 投票:0回答:1

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

我相信这是一个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是53,正如数学所预测的)。这些拟合应该是相同的,但它们不是!我给出了更多的意见。我在下面给出了更多的意见。请你告诉我发生了什么,问题可能出在哪里?

绘制a=1和sigma=sigma_fit的图。

图中sigma=0.65*sigma_fit。


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

请注意一些奇怪的东西......如果我设置sigma = 3* sigma_fit,我得到a = 13,而且拟合度很差(低估)。如果我把它设置为0.2*sigma_fit,我得到的也是一个不好的拟合和a = 10.2 = 5(高估)。以此类推。为什么呢?(顺便说一下,a=1sigma,所以拟合过程应该是有效的)。

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()

然而,如果我手动绘制它,它完全吻合!

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()

EDIT -- 解决了:这是一个绘图问题,需要使用的是

y_fit = piecewise_half_gauss(f_fit, a, sigma = 0.6*sigma_fit)
python numpy matplotlib curve-fitting scipy-optimize
1个回答
0
投票

问题出在绘制和拟合参数上 -- 如果我用不同的参数来拟合 sigma我还需要在绘制部分改变它,当我生成了 y_fit:

# 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)
© www.soinside.com 2019 - 2024. All rights reserved.