python - 在 Python 中为线性拟合(具有两个参数)生成置信区间等高线图

Question

我想根据对任意数据集的最小二乘线性拟合，在 Python 中生成置信区间等高线图。我将 polyfit 函数应用于线性拟合(即 y = mx + c)，由 x、y、yerr 数组上的误差加权，并获得最小卡方值及其对应的线性拟合系数。

从这一点来看，我不知道如何绘制与最佳系数值有 1 西格玛偏差的椭圆。我想在 x 轴上绘制 c，在 y 轴上绘制 m，以及单个 1 sigma 等高线。我一直认为我需要找到卡方函数的反函数(在代码中明确定义)，但这在逻辑上没有意义。

最终，我需要一个 chi^2(m, c) = chi^2_min + 1 形式的椭圆。知道我需要使用什么工具吗？

import numpy as np
import matplotlib.pyplot as plt

# set of x,y values (with y errors) to which a linear fit will be applied
x = np.array([1, 2, 3, 4, 5])
y = np.array([1.7, 2.1, 3.5, 3.2, 4.4])
erry = np.array([0.2, 0.2, 0.2, 0.3, 0.3])

# apply fit to x,y array weighted by 1/erry^2
p2, V = np.polyfit(x, y, 1, w=1/erry, cov=True)

# define a chi square function into which parameter estimates are passed
def chisq(param1, param0):
    csq = np.sum(((param1*x + param0 - y)/erry) ** 2)
    return csq

# arrange labels for the coefficients so matches form y = theta1*x + theta0
theta1 = p2[0]
theta0 = p2[1]
# show coeffs with corresponding stat errors
print("a1 = ",theta1,"+-",np.sqrt(V[0][0]))
print("a0 = ",theta0,"+-",np.sqrt(V[1][1]))

# define arrays for the parameters running between (arbitrarily) parameter +/- 0.3
run1 = np.array([theta1-0.3, theta1-0.2, theta1-0.1, theta1, theta1+0.1, theta1+0.2, theta1+0.3])
run0 = np.array([theta0-0.3, theta0-0.2, theta0-0.1, theta0, theta0+0.1, theta0+0.2, theta0+0.3])

# define the minimum chi square value readily
chisqmin = chisq(run1[4],run0[4])

# Would like to produce a contour at one sigma from min chi square value,
# i.e. obeys ellipse eqn. chi^2(theta0, theta1) = chisqmin + 1

# add lines one sigma away from the optimal parameter values that yield the min chi square value
plt.axvline(x=theta0+np.sqrt(V[1][1]),color='k',linestyle='--')
plt.axvline(x=theta0-np.sqrt(V[1][1]),color='k',linestyle='--')
plt.axhline(y=theta1+np.sqrt(V[0][0]),color='k',linestyle='--')
plt.axhline(y=theta1-np.sqrt(V[0][0]),color='k',linestyle='--')
plt.xlabel(r'$\theta_{0}$')
plt.ylabel(r'$\theta_{1}$')

Answer 1

(您可能在 https://stats.stackexchange.com/ 上运气更好；这是一个统计量很大的问题)

据我了解，您想计算出 χ² 如何随最佳拟合线的梯度和截距(m 和 c)变化。这应该可以通过创建一个可能的 m 和 c 值的数组，为每对计算出 χ² 并绘制这个新的二维数组的等高线来实现。

这是一个基于您的代码的快速示例，它使用 np.linspace() 创建可能的 m 和 c 值的数组，并仅绘制结果卡方的等高线 - 您需要编辑它以获得等高线1 西格玛偏差，但希望这是朝着正确方向迈出的一步。

import numpy as np
import matplotlib.pyplot as plt


# define a chi square function into which parameter estimates are passed
def chisq(param1, param0, x, y, erry):
    csq = np.sum(((param1 * x + param0 - y) / erry) ** 2)
    return csq


def main():
    # set of x,y values (with y errors) to which a linear fit will be applied
    x = np.array([1, 2, 3, 4, 5])
    y = np.array([1.7, 2.1, 3.5, 3.2, 4.4])
    erry = np.array([0.2, 0.2, 0.2, 0.3, 0.3])

    # apply fit to x,y array weighted by 1/erry^2
    p2, V = np.polyfit(x, y, 1, w=1 / erry, cov=True)

    # arrange labels for the coefficients so matches form y = theta1*x + theta0
    theta1 = p2[0]
    theta0 = p2[1]
    # show coeffs with corresponding stat errors
    print("a1 = ", theta1, "+-", np.sqrt(V[0][0]))
    print("a0 = ", theta0, "+-", np.sqrt(V[1][1]))

    # define arrays for parameters running between the mean value +- 1 standard deviation
    run1 = np.linspace(theta1 - np.sqrt(V[0][0]), theta1 + np.sqrt(V[0][0]))
    run0 = np.linspace(theta0 - np.sqrt(V[1][1]), theta0 + np.sqrt(V[1][1]))

    # Work out a 2d array of chi square values for each of the possible m and c values
    chi_square_values = np.array(
        [
            [chisq(run1[i], run0[j], x, y, erry) for j in range(len(run0))]
            for i in range(len(run1))
        ]
    )

    plt.contourf(chi_square_values)
    plt.show()

    print(chi_square_values)


if __name__ == "__main__":
    main()