Figure 9.22:

Monte Carlo evaluation of confidence intervals.

Code for Figure 9.22

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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
# [makes] nthmonte.dat
# [rule] copymat
# Monte Carlo confidence-coverage check; mirrors nthmonte.m structure.
from paresto import Paresto
from misc import save_ascii
import numpy as np
from scipy.stats import f as f_dist

# Model
pe = Paresto()
pe.print_level = 0
pe.transcription = 'shooting'
pe.x = ['ca']
pe.p = ['k', 'n']
pe.d = ['m_ca']
pe.ode = lambda t, y, p: [-p.k * y.ca**p.n]
pe.lsq = lambda t, y, p: [y.ca - y.m_ca]

tfinal = 5.0
nts    = 100
tout   = np.linspace(0, tfinal, nts)
pe.tout = tout

pe.initialize()

# True parameters
kac, nac, ca0ac = 0.5, 2.5, 2.0
thetaac = np.array([kac, nac, ca0ac])  # order matches Paresto's est.thetavec: [k, n, ca]

# Noise-free trajectory
y_ac, _ = pe.simulate({'k': kac, 'n': nac, 'ca': ca0ac})
y_ac = np.array(y_ac).flatten()

measvar    = 1e-2
measstddev = np.sqrt(measvar)
np.random.seed(0)

# Initial guess and bounds
small = 1e-3
large = 5.0
pe.ic.k  = kac;  pe.ic.n = nac;  pe.ic.ca = ca0ac
pe.lb_ub()
pe.lb.k = small;  pe.ub.k = large
pe.lb.n = small;  pe.ub.n = large
pe.lb.ca = small; pe.ub.ca = large

np_     = 3
nmonte  = 10  # change to 500 for figure in text

alph = np.full(nmonte, np.nan)

for i in range(nmonte):
    noise   = measstddev * np.random.randn(nts)
    y_noisy = y_ac + noise

    est, yy, pp = pe.optimize(y_noisy.reshape(1, nts, 1))
    conf = pe.confidence(est, 0.95)

    samplevar = float(est.f) / (nts - np_)
    d = np.array(est.thetavec).flatten() - thetaac
    cont = float(d @ conf.H @ d / (2 * np_ * samplevar))
    alph[i] = float(f_dist.cdf(cont, np_, nts - np_))
    if (i + 1) % 50 == 0:
        print(f'# Monte Carlo iter {i+1}/{nmonte}')

# Coverage curve: actual fraction within each alpha-level ellipse
amin, amax, nas = 0.0, 0.99995, 100
alphavec = np.linspace(amin, amax, nas)
expected = alphavec * nmonte
actual   = np.array([np.sum(alph <= a) for a in alphavec], dtype=float)

table = np.column_stack([alphavec, expected, actual])
save_ascii('nthmonte.dat', table)