Figure 2.3:

Estimated reaction rates from 500 production-rate measurements subject to measurement noise.

Figure 2.3

Code for Figure 2.3

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
# Converted from rate_estimation_lots.m
import numpy as np

stoi = np.array([
    [0, 1, 0, -1, -1, 1],
    [-1, 1, 1, -1, 0, 0],
    [1, 0, -1, 0, -1, 1]
], dtype=float)

r = np.array([1., 2.])
A = stoi[:2, :]
R = A.T @ r

rng = np.random.RandomState(0)

npoints = 500
Rmeas = np.zeros((6, npoints))
for i in range(npoints):
    Rmeas[:, i] = 0.05 * rng.randn(6) + R

rest = (np.linalg.inv(A @ A.T) @ A @ Rmeas).T

with open('rate_estimation_lots.dat', 'w') as myfile:
    for row in rest:
        myfile.write(f'{row[0]:8.4f} {row[1]:8.4f} \n')