Figure 9.13:

Parameter estimates with only 10 data points.

Code for Figure 9.13

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
# Converted from parestonedata.m - single Arrhenius dataset with fit
import numpy as np
from scipy.stats import chi2, f as f_dist
from misc import octave_save

alpha   = 0.95
p       = 2
ndata   = 10
chisq   = chi2.ppf(alpha, p)
Fstat   = f_dist.ppf(alpha, p, ndata - p)
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
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)
data_out   = np.column_stack([1.0/Tmeas, lnkmeas])

nplot  = 20
Tplot  = np.linspace(Tmin, Tmax, nplot)
Xplot  = np.column_stack([np.ones(nplot), -1.0/Tplot + 1.0/np.mean(Tmeas)])
lnkfit = Xplot @ thetacenter
fit    = np.column_stack([1.0/Tplot, lnkfit])

octave_save('parestonedata.dat', ('fit', fit), ('data', data_out))