如何使用 python 在符号图中绘制 2 sigma 变化

Question

我正在尝试使用来自 scipy 的 curve_fit 函数将一些样本数据拟合到一个符号图中。我的最佳拟合曲线在我遵循的代码中看起来不错，但是我遇到了 2 条 sigma 曲线的问题，我想同时显示最佳拟合曲线和灰色填充曲线。我的代码如下所示：

import sys
import os
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import scipy.optimize as optimization

M = np.array([-2, -1, 0, 1, 2, 3,4])
Y_z = np.array([0.05, 0.2, 3, 8, 50, 344, 2400 ])

# curve fit linear function
def line(x, a, b):
    return a*x+b

popt, pcov = curve_fit(line, M, np.log10(Y_z))     # change here

# plotting
plt.semilogy(M , Y_z, 'o')
plt.semilogy(M, 10**line(M, popt[0], popt[1]), ':', label = 'curve-fit')

# plot 1 sigma -error
y1 = 10**(line(M, popt[0] + pcov[0,0]**0.5, popt[1] - pcov[1,1]**0.5))
y2 = 10**(line(M, popt[0] - pcov[0,0]**0.5, popt[1] + pcov[1,1]**0.5))
plt.semilogy(M, y1, ':')
plt.semilogy(M, y2, ':')
plt.fill_between(M, y1, y2, facecolor="gray", alpha=0.15)

plt.xlabel(r"$\log X$")
plt.ylabel('Y')   
plt.legend()
plt.show()

非常感谢您对方差曲线的帮助

Answer 1

原则上，线性拟合根本不需要非线性最小二乘曲线拟合：线性回归应该可以工作。

也就是说，为了解决您的问题，您可能会发现 lmfit (http://lmfit.github.io/lmfit-py/) 在这里很有用。它具有更高级别和更 Pythonic 的曲线拟合方法，并添加了许多功能。其中之一是针对选定的 sigma 值计算结果中的不确定性。

要使用 lmfit，它看起来像

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as optimization

import lmfit

M = np.array([-2, -1, 0, 1, 2, 3,4])
Y_z = np.array([0.05, 0.2, 3, 8, 50, 344, 2400 ])

# curve fit linear function
def line(x, a, b):
    return a*x+b

# set up model and create parameters from model function
# note that function argument names are used for parameters
model = lmfit.Model(line)
params = model.make_params(a=1, b=0)

result = model.fit(np.log10(Y_z), params, x=M)

print(result.fit_report())

这将打印出这样的适合度报告：

[[Model]]
    Model(line)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 8
    # data points      = 7
    # variables        = 2
    chi-square         = 0.10468256
    reduced chi-square = 0.02093651
    Akaike info crit   = -25.4191304
    Bayesian info crit = -25.5273101
[[Variables]]
    a:  0.77630819 +/- 0.02734470 (3.52%) (init = 1)
    b:  0.22311337 +/- 0.06114460 (27.41%) (init = 0)
[[Correlations]] (unreported correlations are < 0.100)
    C(a, b) = -0.447

您可以将最佳拟合结果中的 2-sigma 不确定性计算为

# calculate 2-sigma uncertainty in result
del2 = result.eval_uncertainty(sigma=2, x=M)

然后使用这个和拟合结果来绘制结果（从你的表格稍微修改）：

plt.plot(M, np.log10(Y_z), 'o', label='data')
plt.plot(M, result.best_fit, ':',  label = 'curve-fit')  
plt.fill_between(M, result.best_fit-del2, result.best_fit+del2, facecolor="grey", alpha=0.15)

plt.xlabel(r"$\log X$")
plt.ylabel('Y')
plt.legend()
plt.show()

应该会产生类似的情节

希望有帮助。