Page 430 - Numerical Methods for Chemical Engineering

Bayesian multiresponse regression 419

x 0 = [0.1; 0.1; 0]; % set initial condition
% perform dynamic simulation
[t traj,x traj] = ode45(@batch kinetics dynamics ex, . . .
t span, x 0, [], k1);
% extract response predictions
Y hat=x traj;
return;
batch kinetics dynamics ex.m is the same routine as used previously in the least-squares
fitting of k 1 from the single-response data set c C vs. t.

MCMC simulation with the multiresponse marginal posterior density

−N/2
With the multiresponse marginal posterior density p(θ|Y) ∝|S(θ)| , we use MCMC
simulation to estimate posterior predictions of the form
'
E[g|Y] = g(θ)p(θ|Y)dθ (8.200)
P

We generate, with the Metropolis algorithm, a sequence of states {θ [m] } drawn from p(θ|Y),
and approximate the expectation by

N s
1 [m]
E[g|Y] ≈ g θ (8.201)
N s
m=1
Such calculations are performed by the routine
g pred = Bayes MCMC pred MR(X pred, Y, . . .
fun yhat, fun g, theta 0, MCOPTS, Param);
X pred, Y, fun yhat, and theta 0 take their previous definitions. fun g is the name of a
function that computes g(θ):
g = feval(fun g, theta, Param);

Param is an optional structure that can be passed to fun g.
MCOPTS controls the performance of the MCMC simulation and has the following fields
(and default values):
MCOPTS.N equil: the number of equilibration MC iterations to be performed before sam-
pling begins (1000).
MCOPTS.N samples: the number of MC iterations run during sampling. The larger this
number, the more accurate the prediction (10 000).
MCOPTS.delta theta: the size of the displacements in each parameter during an MC move
relative to average magnitudes of each parameter (0.1). This value is changed to meet
the target fraction of accepted moves.
MCOPTS.accept tar: target fraction of proposed MC moves that are accepted (0.3).