SAS Proc MCMC贝叶斯模型参数约束指定方法问询

问题描述 #

需要估计包含潜变量和两组对比的贝叶斯logistic模型，要求ItemOne1与ItemOne2参数完全相等。当前SAS代码可运行，但结果与频率派方法不符，且ItemOne1、ItemOne2的估计值差异明显。模型中两组theta的先验设定不同：组1的theta均值固定为0，组2的theta均值为带超先验N(0,1)的mu_g。此前尝试用beginnodata/endnodata语法施加约束，但不确定是否正确。

原代码如下：

proc mcmc data=lsat_g outpost=lsat_bayes_1p_unconstrained seed=23 nthreads=1 nbi=5000 nmc=20000;
   array ItemOne[2] ItemOne1 ItemOne2;
   array b1[4]; array b2[4]; array x[5];array mu[2] 0 mu_g;   
   parms ItemOne: 0; parms a1 a2 1; parms b1: b2: 0; 

   prior b1: b2: ~ normal(0, var=16);
   prior a1 a2 ~ lognormal(0, var=9); 

   parms mu_g 0;
   prior mu_g ~ normal(0, var=1);

   beginnodata;
/*      if (ItemOne1 - ItemOne2 = 0.) then lp = lpdfnorm(ItemOne1,0,4);*/
/*  else lp = .;*/
    lp = lpdfnorm(ItemOne1,0,4);
    prior ItemOne1 ItemOne2 ~ general(lp);
   endnodata;

   random theta ~ normal(mu[group], var=1) subject=_obs_; 
   llike=0;
   do j=1 to 5;                     
      if group=1 and j=1 then do;
         prob = logistic(a1*(theta-ItemOne1));
         llike = llike + x[j] * log(prob) + (1 - x[j]) * log(1 - prob);
      end;
      else if group=1 and j > 1 then do;
         prob = logistic(a1*(theta-b1[j-1]));
         llike = llike + x[j] * log(prob) + (1 - x[j]) * log(1 - prob);
      end;
      else if group =2 and j = 1 then do;
                prob = logistic(a2*(theta-ItemOne2));
                llike = llike + x[j] * log(prob) + (1 - x[j]) * log(1 - prob);
      end;
      else if group =2 and j > 1 then do;
                prob = logistic(a2*(theta-b2[j-1]));
                llike = llike + x[j] * log(prob) + (1 - x[j]) * log(1 - prob);
      end;
       
   end;    
   model general(llike);   
run;

问题根源 #

原代码中的约束写法无效：用prior ItemOne1 ItemOne2 ~ general(lp)只是给两个参数指定了相同的先验密度，但参数本身仍会独立采样，并没有强制两者取值相等，这是导致估计值差异大的核心原因。

正确的约束方法 #

方法1：合并参数（推荐，高效稳定） #

直接定义单个参数（如ItemOne），让两组的第一个项目都调用该参数，从根源上保证参数相等，同时减少模型维度，提升采样效率。

修改后的完整代码：

proc mcmc data=lsat_g outpost=lsat_bayes_1p_constrained seed=23 nthreads=1 nbi=5000 nmc=20000;
   array b1[4]; array b2[4]; array x[5];array mu[2] 0 mu_g;   
   parms ItemOne 0; /* 合并为单个参数，强制两组ItemOne相等 */
   parms a1 a2 1; parms b1: b2: 0; 

   prior b1: b2: ~ normal(0, var=16);
   prior a1 a2 ~ lognormal(0, var=9); 
   prior ItemOne ~ normal(0, var=4); /* 给合并后的参数指定先验 */

   parms mu_g 0;
   prior mu_g ~ normal(0, var=1);

   random theta ~ normal(mu[group], var=1) subject=_obs_; 
   llike=0;
   do j=1 to 5;                     
      if group=1 and j=1 then do;
         prob = logistic(a1*(theta-ItemOne)); /* 组1调用ItemOne */
         llike = llike + x[j] * log(prob) + (1 - x[j]) * log(1 - prob);
      end;
      else if group=1 and j > 1 then do;
         prob = logistic(a1*(theta-b1[j-1]));
         llike = llike + x[j] * log(prob) + (1 - x[j]) * log(1 - prob);
      end;
      else if group =2 and j = 1 then do;
                prob = logistic(a2*(theta-ItemOne)); /* 组2同样调用ItemOne */
                llike = llike + x[j] * log(prob) + (1 - x[j]) * log(1 - prob);
      end;
      else if group =2 and j > 1 then do;
                prob = logistic(a2*(theta-b2[j-1]));
                llike = llike + x[j] * log(prob) + (1 - x[j]) * log(1 - prob);
      end;
   end;    
   model general(llike);   
run;

方法2：使用constraint语句（保留原参数名） #

如果需要保留ItemOne1和ItemOne2的参数名，可在beginnodata块中添加constraint语句强制两者相等，同时只需给其中一个参数指定先验即可。

完整代码：

proc mcmc data=lsat_g outpost=lsat_bayes_1p_constrained seed=23 nthreads=1 nbi=5000 nmc=20000;
   array ItemOne[2] ItemOne1 ItemOne2;
   array b1[4]; array b2[4]; array x[5];array mu[2] 0 mu_g;   
   parms ItemOne: 0; parms a1 a2 1; parms b1: b2: 0; 

   prior b1: b2: ~ normal(0, var=16);
   prior a1 a2 ~ lognormal(0, var=9); 
   prior ItemOne1 ~ normal(0, var=4); /* 仅给ItemOne1指定先验 */

   parms mu_g 0;
   prior mu_g ~ normal(0, var=1);

   beginnodata;
      constraint ItemOne1 = ItemOne2; /* 强制两参数相等 */
   endnodata;

   random theta ~ normal(mu[group], var=1) subject=_obs_; 
   llike=0;
   do j=1 to 5;                     
      if group=1 and j=1 then do;
         prob = logistic(a1*(theta-ItemOne1));
         llike = llike + x[j] * log(prob) + (1 - x[j]) * log(1 - prob);
      end;
      else if group=1 and j > 1 then do;
         prob = logistic(a1*(theta-b1[j-1]));
         llike = llike + x[j] * log(prob) + (1 - x[j]) * log(1 - prob);
      end;
      else if group =2 and j = 1 then do;
                prob = logistic(a2*(theta-ItemOne2));
                llike = llike + x[j] * log(prob) + (1 - x[j]) * log(1 - prob);
      end;
      else if group =2 and j > 1 then do;
                prob = logistic(a2*(theta-b2[j-1]));
                llike = llike + x[j] * log(prob) + (1 - x[j]) * log(1 - prob);
      end;
   end;    
   model general(llike);   
run;

额外说明 #

• 若需要对其他b参数（如b1[j]与b2[j]）也施加相等约束，可采用同样逻辑：要么合并参数，要么添加constraint b1[j] = b2[j];语句。
• 方法1减少了模型参数数量，采样速度更快、结果更稳定，优先推荐使用。

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