这需要一些蛮力，但是下面的步骤会在正确的间隔内生成一个椭圆。这不是一个令人满意的解决方案，但它产生了我所追求的情节。在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)
转播
