Figure 9.14:

Confidence intervals with known (solid line) and unknown (dashed line) error variance.

Code for Figure 9.14

Text of the GNU GPL.

main.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
# Converted from parestellipse.m - chi-squared vs F confidence ellipses
import numpy as np
from misc import octave_save, ellipse

p      = 2
ndata  = 10
chisq  = 5.99   # chi^2 95% level, 2 params
Fstat  = 4.46   # F 95% level, 2 params, 8 df
lnk0   = 1.0;  E = 100.0
Tmin   = 300.0; Tmax = 500.0
Tmeas  = np.linspace(Tmin, Tmax, ndata)
X      = np.column_stack([np.ones(ndata), -1.0/Tmeas])
lnk    = X @ np.array([lnk0, E])
measvar = 1e-3
# Match seed used in parestonedata.py / onedata.py / untransform.py so the
# same noise sample drives both the data plot and this ellipse plot
# (the variance must be estimated from the same data shown in the figure).
np.random.seed(382)
lnkmeas = lnk + np.sqrt(measvar) * np.random.randn(ndata)

Tcenter  = -1.0/Tmeas + 1.0/np.mean(Tmeas)
Xcenter  = np.column_stack([np.ones(ndata), Tcenter])
thetacenter = np.linalg.solve(Xcenter.T @ Xcenter, Xcenter.T @ lnkmeas)
thetacentertr = thetacenter.reshape(1, -1)

residual  = lnkmeas - Xcenter @ thetacenter
samplevar = residual @ residual / (ndata - p)

npts  = 181
amat  = Xcenter.T @ Xcenter
lnkest = thetacenter[0]; Eest = thetacenter[1]

# ellipse 1: chi-squared
level1 = chisq * measvar
x1, y1, major1, minor1, bbox1 = ellipse(amat, level1, npts)
x1 += lnkest; y1 += Eest
bbox1[:, 0] += lnkest; bbox1[:, 1] += Eest
outline1 = np.column_stack([x1, y1])

# ellipse 2: F-stat
level2 = Fstat * p * samplevar
x2, y2, major2, minor2, bbox2 = ellipse(amat, level2, npts)
x2 += lnkest; y2 += Eest
bbox2[:, 0] += lnkest; bbox2[:, 1] += Eest
outline2 = np.column_stack([x2, y2])

octave_save('parestellipse.dat',
            ('thetacentertr', thetacentertr),
            ('bbox1', bbox1), ('outline1', outline1),
            ('bbox2', bbox2), ('outline2', outline2))