1、the monte carlo simulation code in SAS: Example1: https://support.sas.com/rnd/app/stat/examples/BayesStd/new_example/index.html

proc mcmc data=class outpost=classout nmc=10000 thin=2 seed=246810;
   parms beta0 0 beta1 0;
   parms sigma2 1;
   prior beta0 beta1 ~ normal(mean = 0, var = 1e6);
   prior sigma2 ~ igamma(shape = 3/10, scale = 10/3);
   mu = beta0 + beta1*height;
   model weight ~ n(mu, var = sigma2);

What does the "parm sigma2 1"mean? Does it mean the sigma2's parameter begins from 1 in the simulation and why it begins from 1 but not 0.1、0.011、0.0001 or others?

Just like Example2 below: https://www.phusewiki.org/docs/Conference%202017%20AS%20Papers/AS03.pdf

the variable "a_Sex" begins from 0.1 and the sigma2 begins from 0.0011. So how to determine the variable's initial value?

proc mcmc data=Orthodont_tr outpost=outpost nbi=5000 nmc=25000 thin=15 plots=ALL; 
array mean[4]; 
array a_age[3]      a_age10 a_age12 a_age14; 
array a_sex_age[3]  a_sex_age10  a_sex_age12  a_sex_age14; 
array dist[4]       dist1 dist2 dist3 dist4; 
parms intercept 0.1 
      a_Sex 0.1 
      a_age: 0 
      a_sex_age: 0; 
parms S2 0.0011; 
parms Rho; 

prior intercept   ~ normal(0,var=1e6); 
prior a_Sex       ~ normal(0,var=1e6); 
prior a_age:      ~ normal(0,var=1e6); 
prior  a_sex_age: ~ normal(0,var=1e6); 
prior S2          ~ igamma(shape=0.01, scale=0.01); 
prior Rho         ~ uniform(-1,1); 

mean[1] = intercept + a_Sex * (Sex=1)                            ;/*time = 8*/ 
mean[2] = intercept + a_Sex * (Sex=1)+ a_age10 +  a_sex_age10 * (Sex=1);/*time = 10*/ 
mean[3] = intercept + a_Sex * (Sex=1)+ a_age12 +  a_sex_age12 * (Sex=1);/*time = 12*/ 
mean[4] = intercept + a_Sex * (Sex=1)+ a_age14 +  a_sex_age14 * (Sex=1);/*time = 14*/ 

model dist ~ mvnar(mean, var=S2, rho=Rho); 


2、Why the sigma2's prior distribution in the 1st example is like this:

prior sigma2 ~ igamma(shape = 3/10, scale = 10/3)

But in the 2rd example:

prior S2 ~ igamma(shape=0.01, scale=0.01)

Why are they different? How to determine the sigma2's priori distribution and why the sigma2 has an igamma distribution?

@merv Do you mean I need to use this formula below to calculate the igamma distribution's "shape(α)" and "scale(1/β)"?

But obviously it is impossible to get the "α=0.01" in Example2 through the "α+n/2"…… enter image description here

  • $\begingroup$ Please don't post image of your code, paste it and use code formatting instead. $\endgroup$
    – Stuart.Sklinar
    Commented Dec 27, 2019 at 15:21
  • $\begingroup$ @Stuart.Sklinar I have pasted all of the code in the Example1 and two key codes about the question in the Example 2. $\endgroup$
    – tumidou
    Commented Dec 27, 2019 at 15:30
  • $\begingroup$ PhUse 2017 are coded images - equally, you may be better off at math.stackexchange.com $\endgroup$
    – Stuart.Sklinar
    Commented Dec 27, 2019 at 15:31
  • 2
    $\begingroup$ @Stuart.Sklinar Ok,thanks. I have updated it. $\endgroup$
    – tumidou
    Commented Dec 27, 2019 at 15:50
  • $\begingroup$ Have you tried changing the parms values? In MCMC sampling for most models, it shouldn't matter where you initialize - as long as it eventually samples from the typical set, it should be good. If it is sensitive to these init values, then that could indicate a problem, or at least a need to discard a larger burn-in set (though the current nbi=5000 is likely sufficient). $\endgroup$
    – merv
    Commented Dec 27, 2019 at 15:50


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.