Page 417 - Numerical Methods for Chemical Engineering

406 8 Bayesian statistics and parameter estimation

MCOPTS.N equil, the number of equilibration MC iterations to be performed before sam-
pling begins (1 000);
MCOPTS.N samples, the number of samples taken to compute the prediction (10,000);
MCOPTS.frac theta, the fraction of Monte Carlo moves that displace θ; the remainder
displace σ (0.5);
MCOPTS.delta theta, the initial relative size of the random Gaussian displacement in each
θ j during proposed move (0.1);
MCOPTS.delta sigma, the initial size of the random Gaussian displacement of
[0]
σ(min{0.1|σ |, √ eps});
MCOPTS.accept tar, the target fraction of moves that are accepted (0.3). The sizes of the
proposed MC moves are adjusted dynamically to achieve this target.

Example. Protein expression data for bacterial strains

We consider again the protein expression data (8.35), where we have fitted the linear model

0, wild-type
y = θ 1 + θ 2 x + ε x = (8.162)
1, mutant
and have obtained the least-squares estimates
θ LS,1 = 113.4 θ LS,2 = 6.7750 (8.163)

We now compute the probability that θ 2 ≥ 5; that is, that the protein expression of the
mutant strain is five units above that of the wild-type strain. We use MCMC simulation to
compute the expectation of the indicator

1, θ 2 ≥ 5
I θ 2 ≥5 = (8.164)
0, θ 2 < 5
We perform this calculation with the code

X pred=[10;10;10;10;11;11;11;11];
y = [121.9; 113.4; 112.2; 106.1; 120.7; 119.5; 116.5; 124.0];
theta 0 = [113.4; 6.7750]; sigma 0=5;
MCOPTS.N equil = 2000; MCOPTS.N samples = 50000;
fun yhat = ‘calc yhat linear model’;
fun g = ‘calc g 1Dinterval’;
Param.j = 2; Param.val lo = 5; Param.val hi = 1e6;
g pred = Bayes MCMC pred SR(X pred, y, fun yhat, fun g,...
theta 0, sigma 0, MCOPTS, Param),
where we have supplied the routines,

% calc yhat linear model.m
function y hat = calc yhat linear model(theta,X);
y hat = X*theta;
return;