Page 416 - Numerical Methods for Chemical Engineering

MCMC computation of posterior predictions 405

generate a new trial state by displacing at random either θ or σ. For some specified fraction
of θ moves f θ ∈ (0, 1), if a random number u, drawn from a uniform distribution on [0, 1],
is less that or equal to f θ , we propose a random displacement of θ:

(new) [m]
θ = θ
j j + θ M j γ j
(8.159)
{γ j } are drawn at random from N(µ = 0,σ = 1)
√
where M j = min{ |θ j | , eps}, |θ j | being the running average of the magnitude of θ j .
Else, we propose a random displacement of σ:

σ (new) = σ [m] + σ γ
(8.160)

γ is drawn at random from N(µ = 0,σ = 1)
θ and σ are relative step sizes that are adjusted throughout the simulation to meet a target
fraction of moves that are accepted (e.g. 0.3).
The probability of accepting the proposed move is
(new) (
π s θ ,σ (new) y
α([m] → (new)) = min 1, (8.161)
π s θ [m] ,σ [m] y

To decide whether to accept the move or not, we generate a random number u drawn
from a uniform distribution on [0, 1]. If u ≤ α([m] → (new)), we accept the move

and (θ [m+1] ,σ [m+1] ) = (θ (new) ,σ (new) ). Else, we reject the move and reuse the old value
(θ [m+1] ,σ [m+1] ) = (θ [m] ,σ [m] ). It is common to run the MCMC simulation for some num-
ber of steps to “equilibrate” the sequence before beginning the sampling stage to compute
the expectation.
The routine Bayes MCMC pred SR.m implements this method to compute the posterior
prediction of a function g(θ,σ) using single-response data, and is invoked by the command
g pred = Bayes MCMC pred SR(X pred, y, fun yhat, fun g, ...
theta 0, sigma 0, MCOPTS, Param);
X pred is a matrix that contains in each row the predictor values for the corresponding
experiment. y is the vector of measured responses. fun yhat is a user-supplied function that
returns the vector of predicted responses:
y hat = feval(fun yhat, theta, X pred);
fun g is a user-supplied function that returns the value of g(θ,σ):

g = feval(fun g, theta, sigma, Param);
Param is an optional structure of parameters that can be passed through
Bayes MCMC pred SR to fun g. theta 0 and sigma 0 are the initial values of θ and σ
that initiate the Markov chain.
MCOPTS is a structure that controls the operation of the MCMC simulation and contains
the following fields (and default values):