Page 435 - Numerical Methods for Chemical Engineering

424 8 Bayesian statistics and parameter estimation

MRSLData.Y 2(:,2) = cB;
MRSLData.Y 2(:,3) = cC;
% specify name of routine that predicts the response
% data matrix for this set of experiments
MRSLData.fun yhat 2 = ‘calc Yhat reaction ex dynamic MR’;

calc Yhat reaction ex dynamic MR.m is the same routine as used previously in fitting the
rate constant using the data of Table 8.3 alone. The following routine predicts the results in
the experiments of Table 8.1 for a trial value of the rate constant:

% calc yhat kinetic ex table1.m
function y hat = calc yhat kinetic ex table1(k 1, X pred);
conc A=X pred(:,1); % [A] in each experiment
conc B=X pred(:,2); % [B] in each experiment
y hat=k 1.*conc A.*conc B;
return;
The cost function and marginal posterior density for the composite data set are computed
by the routine

[F cost, det S, posterior density] = . . .
calc MRSL posterior(theta, MRSLData, det S ref);

For the input value of theta, and the composite data set MRSLData, F cost is the cost function
value, det S is a vector of the |S d (θ)| values in each data set. det S ref is an optional vector
of values |S d (θ (ref) )| used to rescale the posterior density, as above for θ (ref) = θ M . If this
argument is not included, default values of 1 are used.
The following code computes the cost function and |S d (θ)| values for an initial rate
[0]
constant guess of k = 0.001:
1
theta 0 = 0.001;
[F cost 0, det S 0] = calc MRSL posterior(theta 0, MRSLData),
The results are

F cost 0= −203.6059
det S 0 = 1.0e-008 *
0.6012
0.0024
From an initial guess of the parameters, simulated annealing with quenching is used to search
for the global minimum (the returned value is at least a local minimum). This calculation
is performed by the routine

[theta, F cost, det S] = sim anneal MRSL( . . .
MRSLData, theta 0, N iter, . . .
freq quench, freq reset, T 0);