Page 428 - Numerical Methods for Chemical Engineering

Bayesian multiresponse regression 417

Numerically computing the parameter estimate
−N/2
From the marginal posterior density p(θ|Y) ∝|S(θ)| , we identify the most probable
estimate θ M as that minimizing |S(θ)|. We compute |S(θ)| by performing a LU decompo-
sition,

P(θ)S(θ) = L(θ)U(θ) (8.197)

As |S(θ)| > 0 for the positive-definite matrix S(θ), we have

L
0
|S(θ)|= |U ll (θ)| (8.198)
l=1
Here, we have used the fact that |L|= 1 and |P|=±1. As f (x) and ln[ f (x)] have the
same minima, we find θ M by minimizing the cost function

L

F cost (θ) = ln|S(θ)|= ln |U ll (θ)| (8.199)
l=1
which varies more slowly with θ than does |S(θ)|.
Because the gradient vector of this cost function is difficult to determine, we use the
nonlinear simplex method that only requires us to compute the cost function value for each
trial θ. This technique returns only a local minimum; therefore, we augment it by simulated
annealing to search for a global minimum, as is done by the routine sim anneal MR.m:
[theta, det S, det S 0] = sim anneal MR( . . .
X pred, Y, fun yhat, theta 0, N iter, . . .
freq quench, freq reset, T 0);

X pred is the N × M matrix that contains in each row the predictors for the corresponding
experiment. Y is an N × L matrix containing the responses for each experiment in the
corresponding row. fun yhat is the name of a user-supplied routine that returns the predicted
response data matrix for input values of theta (the parameters) and X pred:
Y hat = feval(fun yhat, theta, X pred);

theta 0 contains the initial guesses of the parameters. N iter is the total number of iterations
during the simulated annealing run. freq quench is a number in [0, 1] that sets the frequency
with which the simulation is “quenched” at random times during the simulation. Here,
“quenching” refers to running the nonlinear simplex minimizer fminsearch to find a local
minimum of the cost function. freq reset specifies the frequency with which the state is
reset, at random times, to the best one (lowest cost function) found to date. T 0 is the initial
[0]
annealing temperature. If this argument is not included, a default value of 10|F cost (θ )| is
used.
The routine returns the theta value with the lowest cost function found during the
simulation. det S is the value of |S(θ)| for this optimal parameter value (at least a local
[0]
minimum). det S 0 is the value of |S(θ )|.