Simulation using PROC MCMC

[JOADUTWUM]

Hello all,

I want to simulate my data (a 6 by 6 matrix) using PROC MCMC.

The codes are given below:

proc iml;
N=50000;
Mean = {0.673 0.663 0.625 0.651 0.680 0.704};
cov = {0.003368 0.000770 0.000460 0.0001367 0.000609 0.000474,
0.000770 0.003885 0.000716 -0.00021 0.000540 0.000602,
0.000460 0.000716 0.002478 0.001006 0.001013 0.000737,
0.0001367 -0.00021 0.001006 0.005104 0.001282 0.000681,
0.000609 0.000540 0.001013 0.001282 0.002611 0.001241,
0.000474 0.000602 0.000737 0.000681 0.001241 0.003422};
call randseed(1);
x=Randnormal(N,Mean,Cov);
SampleMean=x[:,];
n=nrow(x);
y= x-repeat(SampleMean, n);
SampleCov = y`* y / (n-1);
print SampleMean Mean, SampleCov Cov;
cname={"x1", "x2", "x3", "x4", "x5", "x6"};
create inputdata from x [colname=cname];
append from x;
close inputdata;
quit;

proc corr data=inputdata plots(maxpoints=NONE)=matrix(histogram);
var x:;
run;
proc mcmc data=inputdata seed=17 nmc=50000 diag=none ;
ods select PostSumInt;
array data[6] x1 x2 x3 x4 x5 x6;
array mu[6] (0 0 0 0 0 0);
array Sigma[2,2,2,2,2,2];
array mu0[6] (0 0 0 0 0 0);
array Sigma0 [6,6] (100 0 0 0 0 100);
array S[6,6] (1 0 0 0 0 1);
parm mu Sigma ;
prior mu ~ mvn(mu0, Sigma0);
prior Sigma ~ iwish(6,S);
rho = sigma[1,2]/sqrt(sigma[1,1]*sigma[2,2]);
model data ~ mvn(mu, Sigma);
run;

ERROR: In the inverse Wishart distribution, Sigma must be a 2-dimensional array. A 6-dimensional array was specified.

I would appreciate your help. 

Thannks 

[Rick_SAS · Accepted solution]

In general, if you are using code that you copied from somewhere, please provide a link to the original source. That is helpful in trying to determine what the code is doing.

In this case, it appears that you copied code from the PROC MCMC documentation. The documentation example is for a bivariate distribution, so the covariance parameters are 2 x 2. Your example uses 6 x 6 covariance matrices. Thus, you need to modify the PROC MCMC code to use 6 x 6 matrices.  However, the example runs on a simulated data set that contains N=50,000 observations, which is quite large. In the following statements, I set N=10,000 and reduced NMC=5000 so that it runs faster.

proc iml;
N=10000;
Mean = {0.673 0.663 0.625 0.651 0.680 0.704};
cov = {0.003368 0.000770 0.000460 0.0001367 0.000609 0.000474,
0.000770 0.003885 0.000716 -0.00021 0.000540 0.000602,
0.000460 0.000716 0.002478 0.001006 0.001013 0.000737,
0.0001367 -0.00021 0.001006 0.005104 0.001282 0.000681,
0.000609 0.000540 0.001013 0.001282 0.002611 0.001241,
0.000474 0.000602 0.000737 0.000681 0.001241 0.003422};
call randseed(1);
x=Randnormal(N,Mean,Cov);

cname={"x1", "x2", "x3", "x4", "x5", "x6"};
create inputdata from x [colname=cname];
append from x;
close inputdata;
quit;
/*
proc corr data=inputdata plots(maxpoints=NONE)=matrix(histogram);
var x:;
run;
*/
proc mcmc data=inputdata seed=17 nmc=5000 diag=none ;
*ods select PostSumInt;
array data[6] x1 x2 x3 x4 x5 x6;
array mu[6] (0 0 0 0 0 0);
array Sigma[6,6];
array mu0[6] (0 0 0 0 0 0);
array Sigma0 [6,6] (100  0   0   0   0   0
                    0  100   0   0   0   0
                    0    0 100   0   0   0
                    0    0   0 100   0   0
                    0    0   0   0 100   0
                    0    0   0   0   0 100);
array S[6,6] (1 0 0 0 0 0
              0 1 0 0 0 0
              0 0 1 0 0 0
              0 0 0 1 0 0
              0 0 0 0 1 0
              0 0 0 0 0 1);
parm mu Sigma ;
prior mu ~ mvn(mu0, Sigma0);
prior Sigma ~ iwish(6,S);
rho = sigma[1,2]/sqrt(sigma[1,1]*sigma[2,2]);  /* ?? not used */
model data ~ mvn(mu, Sigma);
run;

[JOADUTWUM]

Thank you @Rick_SAS 

It did work.