我想在Python中生成一个基于对任意数据集的最小二乘线性拟合的置信区间等高线图。我将多元拟合函数应用于由x、y、yerr阵列上的误差加权的线性拟合(即y=mx+c(，并获得最小卡方值及其对应的线性拟合系数。
从这一点来看，我不知道如何绘制与最佳系数值偏差为1西格玛的椭圆。我想在x轴上画c，在y轴上画m，画一个1西格玛的等值线。我一直认为我需要找到卡方函数的逆函数(在代码中明确定义(，但这在逻辑上没有意义。
最后，我需要一个形式为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}$')
https://stats.stackexchange.com/；这是一个很重统计数据的问题(
据我所知，你想弄清楚2随最佳拟合线的梯度和截距(m和c(而变化。这应该可以通过创建一个可能的m和c值的数组来实现，计算出χ2，并绘制这个新的2d阵列的轮廓。
这里有一个基于您的代码的快速示例，它使用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()
这需要一些蛮力，但以下操作会在正确的间隔内生成一个椭圆。这不是一个令人满意的解决方案，但它产生了我想要的情节。
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)