请教SAS问题 PROC MCMC

rossinifa

本人想用SAS进行Bayesian Logistic-Normal Model

在帮助文件中，（PROC MCMC Example 52.5 Random-Effects Models）
有这样一段程序
proc mcmc data=seeds outpost=postout seed=332786 nmc=100000 thin=10
          ntu=3000 monitor=(beta0-beta3 v);
   ods select PostSummaries ess;
   array delta[21];

   parms delta: 0;
   parms beta0 0 beta1 0 beta2 0 beta3 0 ;
   parms v 1;

   beginnodata;
   sigma = sqrt(v);
   endnodata;
   w = beta0 + beta1*seed + beta2*extract + beta3*seed*extract;
   if ind eq 1 then
      lp = lpdfnorm(delta[ind], w, sigma);
   else
      lp = lp + lpdfnorm(delta[ind], w, sigma);

   prior v   ~ general(-log(v));
   prior beta: ~ general(0);
   prior delta: ~ general(lp);

   pi = logistic(delta[ind]);
   model r ~ binomial(n = n, p = pi);
run;
我个人觉得  prior delta: ~ general(lp);  是不准确 的。
在模型中的random term 如果看成是普通的参数，但是赋予先验分布，那么，理论上应该和把它当成random effect 是等价 的。
我认为应该可以利用prior 语句这样做：
proc mcmc data=seeds outpost=postout seed=332786 nmc=100000 thin=10                                                                     
          ntu=3000 monitor=(beta0-beta3 v);                                                                                             
                                                                                                                                       
   array delta[21];                                                                                                                     
                                                                                                                                       
   parms delta: 0;                                                                                                                     
   parms beta0 0 beta1 0 beta2 0 beta3 0 ;                                                                                             
   parms v 1;                                                                                                                           
                                                                                                                                       
   beginnodata;                                                                                                                        
   sigma = sqrt(v);                                                                                                                     
   endnodata;                                                                                                                           
                                                                                                                                       
                                                                                                                                       
   prior v   ~ general(-log(v));                                                                                                        
   prior beta: ~ general(0);                                                                                                            
   prior delta: ~ normal(0,var=v);                                                                                                      
                                                                                                                                       
w = beta0 + beta1*seed + beta2*extract + beta3*seed*extract;                                                                           
                                                                                                                                       
   pi = logistic(w+delta[ind]);                                                                                                         
   model r ~ binomial(n = n, p = pi);                                                                                                   
run;                                                                                                                                    
        


但是和帮助文件中程序的结果很不一样。
不知道我是错在哪里了？？

请高手帮忙。