MCMC algorithm for chemical kinetic parameter estimation

Question

I'm currently envolved in application of MCMC algorithm to kinetic parameter estimation concerning the example reported below as part of the following reference.

Beers, Kenneth. Numerical Methods for Chemical Engineering: Applications in MATLAB®. New York, NY: Cambridge University Press, November 2006

MatLab Code at page. 411.

% input predictor and response data 
X_pred = [0.1 0.1; 0.2 0.1; 0.1 0.2; 0.2 0.2; 0.05 0.2; 0.2 0.05];
y = [0.0246e-3; 0.0483e-3; 0.0501e-3; 0.1003e-3; 0.0239e-3; 0.0262e-3];

% provide name of routine that returns predicted responses
fun_yhat = 'calc_yhat_kinetic_ex1';
% provide initial guess of parameters
theta_0 = [0.0025 1.0 1.0];
% compute sample variance with initial guess, and use its
% square root as initial guess for sigma
% y_hat = feval(fun_yhat,theta_0,X_pred);
y_hat = feval(fun_yhat,theta_0,X_pred);
RSS_0 = dot(y-y_hat,y-y_hat);
sample_var_0 = RSS_0/(length(y)-length(theta_0));
sigma_0 = sqrt(sample_var_0);
% select the parameters whose 2-D marginal density
% is desired
i_plot_2D = 2; j_plot_2D = 3;
% set the histogram properties
val_lo = [0.8; 0.8]; val_hi = [1.2; 1.2]; N_bins = 50;

%%
% perform the MCMC simulation to compute the
% 2-D marginal posterior density
MCOPTS.N_equil = 50000; % # of equilibration iterations
MCOPTS.N_samples = 25000000; % # of samples
[bin_2Dic, bin_2Djc, bin_2Dp] = ...
    Bayes_MCMC_2Dmarginal_SR(X_pred,y,...
    fun_yhat,i_plot_2D,j_plot_2D,val_lo,val_hi,...
    N_bins,theta_0,sigma_0,MCOPTS);
% generate from the results the 95% HPD region
alpha = 0.05;
HPD_2D = Bayes_2D_HPD_SR(bin_2Dic, bin_2Djc, bin_2Dp, ...
i_plot_2D, j_plot_2D, alpha);

Now, the pre-defined functions are the following ones:

fun_yhat = ‘calc yhat kinetic ex1’

% calc_yhat_kinetic_ex1.m
function y_hat = calc_yhat_kinetic_ex1(theta, X_pred);

N = size(X_pred,1);  % # of experiments

% extract predictors
conc_A = X_pred(:,1);  % [A] in each experiment
conc_B = X_pred(:,2);  % [B] in each experiment

% extract parameters
k_1 = theta(1);  % rate constant
nu_a = theta(2);  % exponent for [A]
nu_b = theta(3); % exponent for [B]

% make predictions
y_hat = zeros(N,1);
for k=1:N
y_hat(k) = k_1*conc_A(k)^nu_a*conc_B(k)^nu_b;
end
% NOTE, WE COULD SIMPLY USE
% yhat = k1.*(conc_A.^nu_A).*(conc_B.^nu_B);
% TO AVOID SLOW FOR LOOP

return;

Bayes_MCMC_2Dmarginal_SR

function [bin_ic, bin_jc, bin_p] = ...
    Bayes_MCMC_2Dmarginal_SR( X_pred, y, ...
        fun_yhat, i_plot, j_plot, val_lo, val_hi, ...
       N_bins, theta_0, sigma_0, MCOPTS);

% initialize the output structures
bin_ic = zeros(N_bins,1);  bin_jc = bin_ic;
bin_delta = (val_hi-val_lo)./N_bins;
for k=1:N_bins
    bin_ic(k) = val_lo(1) + (k-0.5)*bin_delta(1);
    bin_jc(k) = val_lo(2) + (k-0.5)*bin_delta(2);
end
bin_p = zeros(N_bins,N_bins);

% set name of routine that returns the
% appropriate predictors for use in
% generating a 2-D marginal density
fun_g = 'calc_g_2Dinterval';
Param.val_lo = val_lo;
Param.val_hi = val_hi;
Param.i_plot = i_plot;
Param.j_plot = j_plot;
Param.N_bins = N_bins;
Param.bin_delta = bin_delta;

% Now, perform the Markov Chain Monte Carlo
% (MCMC) simulation
g_pred = Bayes_MCMC_pred_SR(X_pred, y, ...
    fun_yhat, fun_g, theta_0, sigma_0, MCOPTS, Param);

% From the results of the MCMC simulation, generate
% the 2-D marginal density.
area_bin_inv = 1/bin_delta(1)/bin_delta(2);
for ki = 1:N_bins
    for kj = 1:N_bins
        label = (ki-1)*N_bins + kj;
        bin_p(ki,kj) = g_pred(label)*area_bin_inv;
    end
end

% make contour and surface plots of the 2-D marginal density
% contour plot
figure;
contour(bin_ic,bin_jc,bin_p);
xlabel(['\theta(', int2str(i_plot),')']);
ylabel(['\theta(', int2str(j_plot), ')']);
title('2-D marginal density from MCMC simulation');
% filled contour plot
figure;
contourf(bin_ic,bin_jc,bin_p);  colorbar;
xlabel(['\theta(', int2str(i_plot),')']);
ylabel(['\theta(', int2str(j_plot), ')']);
title('2-D marginal density from MCMC simulation');
% surface plot
figure;
surf(bin_ic,bin_jc,bin_p);
xlabel(['\theta(', int2str(i_plot),')']);
ylabel(['\theta(', int2str(j_plot), ')']);
zlabel(['p(\theta(', int2str(i_plot), ', ', ...
        int2str(j_plot), '| y)']);
title('2-D marginal density from MCMC simulation');

return;

Bayes_MCMC_pred_SR

function g_pred = Bayes_MCMC_pred_SR(X_pred, y, ...
    fun_yhat, fun_g, theta_0, sigma_0, MCOPTS, Param);

routine_name = 'Bayes_MCMC_pred_SR';
num_arg_in = nargin;

% extract dimension of problem
P = length(theta_0);
N = length(y);

% add default values to MCOPTS as needed
% .N_equil
try
    val = MCOPTS.N_equil;
    if(val < 0)
        error([routine_name, ...
                ': MCOPTS.N_equil < 0']);
    end
catch
    MCOPTS.N_equil = 1000;
end

%  .N_samples
try
    val = MCOPTS.N_samples;
    if(val < 0)
        error([routine_name, ...
                ': MCOPTS.N_samples < 0']);
    end
catch
    MCOPTS.N_samples = 10000;
end

%  .frac_theta
try
    val = MCOPTS.frac_theta;
    if(or((val<=0),(val>=1)))
        error([routine_name, ...
                ': MCOPTS.frac_theta out of (0,1)']);
    end
catch
    MCOPTS.frac_theta = 0.5;
end
        
% .delta_theta
try
    val = MCOPTS.delta_theta;
catch
    MCOPTS.delta_theta = 0.1;
end

% .delta_sigma
try
    val = MCOPTS.delta_sigma;
catch
    MCOPTS.delta_sigma = 0.1*abs(sigma_0);
    if(MCOPTS.delta_sigma < sqrt(eps))
        MCOPTS.delta_sigma = sqrt(eps);
    end
end
    
% .accept_tar
try
    val = MCOPTS.accept_tar;
    if(or((val<=0),(val>=1)))
        error([routine_name, ...
                ': MCOPTS.accept_tar out of (0,1)']);
    end

catch
    MCOPTS.accept_tar = 0.3;
end

% Now, we initialize the counters that keep track of the
% acceptance fractions and the numbers of samples
% taken.
num_theta_try = 0;
num_theta_accept = 0;
num_sigma_try = 0;
num_sigma_accept = 0;
num_samples_taken = 0;

% initialize the running sums of the magnitudes of the
% parameters, and the values avg_mag used in setting
% the step sizes for each parameter.
sum_mag = abs(theta_0);
avg_mag = sum_mag;
k_small = find(avg_mag < sqrt(eps));
avg_mag(k_small) = sqrt(eps);

% with initial theta_0, sigma_0, compute the
% initial vector of predictions, sum of squared
% errors, and posterior probability.
theta = theta_0;
sigma = sigma_0;
y_hat = feval(fun_yhat,theta,X_pred);
residuals = y - y_hat;
SSE = dot(residuals,residuals);
posterior_density = sigma^(-(N+1))* ...
    exp(-SSE/2/sigma^2);

% allocate space to store prediction of g
% and set to initially contain zeros
if(nargin >= 8)
    g_pred = feval(fun_g,theta,sigma,Param);
else
    g_pred = feval(fun_g,theta,sigma);
end
g_pred = zeros(size(g_pred));

% begin Monte Carlo iterations
delta_theta = MCOPTS.delta_theta;
delta_sigma = MCOPTS.delta_sigma;
N_iter_tot = MCOPTS.N_equil + MCOPTS.N_samples;
for iter_MC = 1 : N_iter_tot

    % occasionally, adjust step sizes to fit
    % target acceptance fraction if we have
    % collected enough data since last
    % adjustment. If the acceptance fraction
    % is too low, decrease the step size of
    % the moves.  If too high, increase the
    % step size of the moves.
    if(num_theta_try > 500)
        frac_accept = ...
            num_theta_accept/num_theta_try;
        if(frac_accept < MCOPTS.accept_tar)
            delta_theta = delta_theta * 0.9;
        elseif(frac_accept > MCOPTS.accept_tar)
            delta_theta = delta_theta * 1.1;
        end
        % reset counters
        num_theta_try = 0;
        num_theta_accept = 0;
    end
    % do the same thing for the sigma moves
    if(num_sigma_try > 500);
        frac_accept = ...
            num_sigma_accept/num_sigma_try;
        if(frac_accept < MCOPTS.accept_tar)
            delta_sigma = delta_sigma * 0.9;
        elseif(frac_accept > MCOPTS.accept_tar)
            delta_sigma = delta_sigma * 1.1;
        end
        num_sigma_try = 0;
        num_sigma_accept = 0;
    end        
    
    % decide whether to do theta move or sigma
    % move
    u = rand;  % generate uniform random # on [0,1]
    if(u <= MCOPTS.frac_theta)  % do theta move
        do_theta_move = 1;
        do_sigma_move = 0;
    else
        do_theta_move = 0;
        do_sigma_move = 1;
    end
    if(do_theta_move)   % do theta move
        theta_new = theta + ...
            delta_theta.*avg_mag.* ...
            randn(size(theta));
        num_theta_try = num_theta_try + 1;
        % compute new predictions with this theta, and
        % the new weight according to the posterior
        % density
        y_hat_new = feval(fun_yhat,theta_new,X_pred);
        residuals_new = y - y_hat_new;
        SSE_new = dot(residuals_new,residuals_new);
        posterior_density_new = sigma^(-(N+1))* ...
            exp(-SSE_new/2/sigma^2);
        % compute the acceptance probability
        accept_prob = exp( (SSE - SSE_new) / (2*sigma^2) );
        if(accept_prob > 1)
            accept_prob = 1;
        end
        % decide whether to accept the move or not
        u = rand;
        if(u <= accept_prob)  % accept the move
            theta = theta_new;
            y_hat = y_hat_new;
            residuals = residuals_new;
            SSE = SSE_new;
            posterior_density = posterior_density_new;
            num_theta_accept = num_theta_accept + 1;
        end
    end
    if(do_sigma_move) % do sigma move
        sigma_new = sigma + delta_sigma*randn;
        num_sigma_try = num_sigma_try + 1;
        % compute the acceptance probability
        if(sigma_new <= 0)
            accept_prob = 0;
        else
            % compute new posterior density (but note that
            % the value of SSE does not change)
            posterior_density_new = -sigma_new^(-(N+1))* ...
                exp(-SSE/2/sigma_new^2);
            accept_prob = (sigma_new/sigma)^(-(N+1)) * ...
                exp(SSE/2*(-1/(sigma_new^2) + 1/(sigma^2)));
            if(accept_prob > 1)
                accept_prob = 1;
            end
        end
        % decide whether to accept the move
        u = rand;
        if(u <= accept_prob)
            sigma = sigma_new;
            posterior_density = posterior_density_new;
            num_sigma_accept = num_sigma_accept + 1;
        end
    end
    
    % if we are in the sampling stage, then measure the function
    % g(theta,sigma) and add to running sum.
    if(iter_MC > MCOPTS.N_equil)
        if(nargin >= 8)
            g_meas = feval(fun_g,theta,sigma,Param);
        else
            g_meas = feval(fun_g,theta,sigma);
        end
        g_pred = g_pred + g_meas;
        num_samples_taken = num_samples_taken + 1;
    end
    
    % increment the values of sum_mag and avg_mag
    sum_mag = sum_mag + abs(theta);
    avg_mag = sum_mag/(iter_MC+1);
    k_small = find(avg_mag < sqrt(eps));
    avg_mag(k_small) = sqrt(eps);
    
end

% get prediction from sum of measured predictions
g_pred = g_pred/num_samples_taken;


return;

I am not able to write the missing function fun_g = 'calc_g_2Dinterval'; enclosed in the Bayes_MCMC_2Dmarginal_SR function. I have tried to write it but I think that some mistake are present because the main code dose not work/run.

function g = calc_g_2Dinterval(theta,sigma,Param);

if((theta(Param.j) >= Param.val_lo(2)) & ...
        (theta(Param.j) <= Param.val_hi(2)))
    g(2) = 1;
else
    g(2) = 0;
end

if((theta(Param.i) >= Param.val_lo(1)) & ...
        (theta(Param.i) <= Param.val_hi(1)))
    g(1) = 1;
else
    g(1) = 0;
end

return;

Do you have any suggestions or solution? Thanks.

The error messages that MatLab prints, are the following one:

Reference to non-existent field 'j'.

Error in calc_g_2Dinterval (line 3)
if((theta(Param.j) >= Param.val_lo(2)) & ...

Error in Bayes_MCMC_pred_SR (line 175)
    g_pred = feval(fun_g,theta,sigma,Param);

Error in Bayes_MCMC_2Dmarginal_SR (line 72)
g_pred = Bayes_MCMC_pred_SR(X_pred, y, ...

Error in example_411 (line 29)
    Bayes_MCMC_2Dmarginal_SR(X_pred,y,...

---

I posted this question on Cross Validated since topic related to statistics, but it was pointed out to me that I had posted on the wrong section. I hope this section (related to programming) is the correct one and I am not off-topic.

I link here the original question: The original question posted on Cross Validated