03-26-2015 02:01 PM
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.
|Header 1||Initial Value||Proc MCMC||Inverse Weibull S. Function|
/* Simulate Weibull data using Proc MCMC */
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);
proc lifereg data=simout_weibull;
/* simulate Weibull using survival function S=exp(-(x/scale)**shape) */
do i=1 to 100000;
proc lifereg data=simout_weibull2;
03-27-2015 11:08 AM
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
03-27-2015 11:34 AM
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.
proc mcmc data=x outpost=simoutweib seed=23 nmc=100000
statistics=(summary interval) ;
ods exclude nobs;
prior x ~ weibull(0,&shape,&scale);
proc lifereg data=simoutweib;