修复MCMC化学动力学参数估计MATLAB缺失函数报错问题

问题背景 #

基于Kenneth Beers所著《Numerical Methods for Chemical Engineering: Applications in MATLAB®》第411页的示例，在应用MCMC算法开展化学动力学参数估计时，编写的calc_g_2Dinterval函数运行出现字段不存在的报错，导致MCMC仿真无法正常执行。

主示例代码 #

% 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);

报错的calc_g_2Dinterval函数 #

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;

错误信息 #

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,...

问题分析与修复方案 #

报错根源是函数中错误引用了不存在的字段Param.j和Param.i，实际Bayes_MCMC_2Dmarginal_SR传递的参数字段为Param.j_plot和Param.i_plot。修复后的函数如下：

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

% 检查第i_plot个参数是否在第一个区间内
if((theta(Param.i_plot) >= Param.val_lo(1)) && ...
        (theta(Param.i_plot) <= Param.val_hi(1)))
    g(1) = 1;
else
    g(1) = 0;
end

% 检查第j_plot个参数是否在第二个区间内
if((theta(Param.j_plot) >= Param.val_lo(2)) && ...
        (theta(Param.j_plot) <= Param.val_hi(2)))
    g(2) = 1;
else
    g(2) = 0;
end

return;

修复说明 #

• 将Param.j替换为Param.j_plot，Param.i替换为Param.i_plot，匹配调用时传递的参数字段
• 调整判断顺序，先处理i_plot对应第一个区间，再处理j_plot对应第二个区间，逻辑更清晰
• 用&&替代&，实现短路逻辑判断，提升效率

内容的提问来源于stack exchange，提问作者andreamav