Figure 9.10:

Distribution of estimated parameters.

Code for Figure 9.10

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
# Converted from arrhenius.m - Arrhenius parameter estimation, confidence ellipse
import numpy as np
from scipy.stats import chi2
from misc import ellipse, octave_save

chisq   = 5.99
lnk0    = 1.;   E = 100.
ndata   = 10;   Tmin = 300.;   Tmax = 500.
Tmeas   = np.linspace(Tmin, Tmax, ndata)
X       = np.column_stack([np.ones(ndata), -1./Tmeas])
lnk     = X @ np.array([lnk0, E])
measvar = 1e-3

np.random.seed(0)
nexpts  = 500
lnkmeas = lnk[:,None] + np.sqrt(measvar)*np.random.randn(ndata, nexpts)

theta = (np.linalg.lstsq(X, lnkmeas, rcond=None)[0]).T  # (nexpts, 2)

Tcenter   = -1./Tmeas + 1./np.mean(Tmeas)
Xcenter   = np.column_stack([np.ones(ndata), Tcenter])
thetacenter = (np.linalg.lstsq(Xcenter, lnkmeas, rcond=None)[0]).T

npts = 181
amat = X.T @ X / measvar
x, y, major, minor, bbox = ellipse(amat, chisq, npts)
x = x + lnk0;  y = y + E
minor[:,0] += lnk0;  minor[:,1] += E
major[:,0] += lnk0;  major[:,1] += E
bbox[:,0]  += lnk0;  bbox[:,1]  += E
elips = np.column_stack([x, y])

octave_save('arrhenius.dat',
            ('theta', theta),
            ('bbox',  bbox),
            ('elips', elips))