在 Python 中为线性拟合（带有两个参数）生成置信区间等高线图

Question

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

从这一点来看，我不知道如何绘制椭圆，使其与最佳系数值有 1 sigma 的偏差。我想在 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 sigma 偏差，但希望这是朝着正确方向迈出的一步。

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

Answer 2

它采取了一些蛮力，但以下会在正确的间隔内产生一个椭圆。这不是一个令人满意的解决方案，但它产生了我所追求的情节。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
import math

# 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])
ax = plt.subplot(111)

# 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 +/- sigma
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]))

# define the minimum chi square value readily
chisqmin = chisq(theta0, theta1)
print(chisqmin)

# 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='--', linewidth=0.8)
plt.axvline(x=theta0-np.sqrt(V[1][1]),color='k',linestyle='--', linewidth=0.8)
plt.axhline(y=theta1+np.sqrt(V[0][0]),color='k',linestyle='--', linewidth=0.8)
plt.axhline(y=theta1-np.sqrt(V[0][0]),color='k',linestyle='--', linewidth=0.8)
plt.plot(theta0, theta1, 'o', markersize=4, color='k')
plt.annotate(r'LS estimate',
            xy=(theta0, theta1), xytext=(-80, -40), textcoords='offset points', fontsize=14,
            arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
plt.annotate(r'$\chi^{2}(\theta_{0}, \theta_{1})$ = $\chi^{2}_{min}$ + 1',
            xy=(1.2, 0.7), xytext=(-22, +30), textcoords='offset points', fontsize=14,
            arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
plt.xlabel(r'$\theta_{0}$', fontsize=16)
plt.ylabel(r'$\theta_{1}$', fontsize=16)
plt.xlim(theta0-2*np.sqrt(V[1][1]), theta0+2*np.sqrt(V[1][1]))
plt.ylim(theta1-2*np.sqrt(V[0][0]), theta1+2*np.sqrt(V[0][0]))

sig0 = np.sqrt(V[1][1])
sig1 = np.sqrt(V[0][0])
rho = V[0][1]/(sig1*sig0)
tantwophi = 2*rho*sig1*sig0/(sig0**2-sig1**2)
twophi = math.atan(tantwophi)
phi = twophi/2
phideg = math.degrees(phi)


ellipse=Ellipse((theta0, theta1), width=2.1*np.sqrt(V[1][1]), 
                height=0.8*np.sqrt(V[0][0]), angle=phideg, color='k', ls='-', lw=1.5, fill=False)
ax.add_patch(ellipse)