Why Aren't The Simulated Data by Proc MCMC Close to The Initial Values?

[SASLEE]

I simulated two data sets of the Weibull distribution with the same initial values. Each data set has 100,000 observations.

The first data set was simulated using Proc MCMC.

The second data set was simulated by inversing survival function of the Weibull (exp(-(x/scale)**shape)).

Then, Proc Lifereg is applied in finding the parameter estimates of the two data sets.

The following is the comparison. Apparently, the data simulated by Proc MCMC are not close to the initial values.

Could you please let me know what option I need to add to Proc MCMC to improve the results?

Thanks you very much.

Lee

Header 1	Initial Value	Proc MCMC	Inverse Weibull S. Function
Scale	10	11.3552	10.0532
Shape	0.5	0.5457	0.5002

%let shape=0.5;
%let scale=10;

/* Simulate Weibull data using Proc MCMC */
data null;
run;
proc mcmc data=null outpost=simout_weibull seed=5678 nmc=100000;
parm x 15;
logpdf_weibull=logpdf('weibull', x, &shape, &scale);
prior x ~ general(logpdf_weibull);
model general(0);
run;

proc lifereg data=simout_weibull;
model x=/dist=weibull;
quit;

/* simulate Weibull using survival function S=exp(-(x/scale)**shape) */
data simout_weibull2;
do i=1 to 100000;
x=&scale*(-log(ranuni(5678)))**(1/&shape);
output;
end;
keep x;
run;

proc lifereg data=simout_weibull2;
model x=/dist=weibull;
quit;

[Rick_SAS]

I'm not an expert in MCMC, so hopefully a real expert will weight in. However, you should realize that the second method draws 100,000 INDEPENDENT observations, whereas the MCMC method draws 100,000 CORRELATED observations. For the MCMC sample, the effective sample size (ESS) is much smaller than 100,000. If you look at the MCMC diagnostics, it looks like the Markov chain has only marginal mixing and that the ESS is less than 1,500.

See the SAS/STAT doc for how to assess Markov chain convergence.

The PROC MCMC documentation has some information about possible ways to handle slow convergence

[lvm]

Here is a better way of doing this, directly using the Weibull distribution function in MCMC. This is based on example 1 in the User's Guide. This gives very good mixing and uncorrelated values. The parameter estimates are very close to the theoretical.

data x;

run;

proc mcmc data=x outpost=simoutweib seed=23 nmc=100000

     statistics=(summary interval) ;

   ods exclude nobs;

   parm x;

   prior x ~ weibull(0,&shape,&scale);

   model general(0);

run;

proc lifereg data=simoutweib;

model x=/dist=weibull;

quit;