* Normal - t data;
data ndata;
infile "N:\Private\Shared\Research\BayesianSEM\GCM with flexible error\Simulation\tdata.txt";
input time id y;
run;

PROC MCMC DATA=ndata NMC=40000  NBI=5000 THIN=2
DIC SEED=17 INIT=RANDOM DIAG=(ESS GEWEKE(F1=.2 F2=.5)) 
STATISTICS(ALPHA=0.05)=(SUMMARY INTERVAL)
OUTPOST=histnorm; 
ODS SELECT PARAMETERS POSTSUMMARIES 
   TADPANEL GEWEKE ESS DIC POSTINTERVALS;
ARRAY b[2] L S;
ARRAY beta[2];
ARRAY Sigma_b[2,2];
ARRAY beta0[2] (0 0);
ARRAY sigma0[2,2] (1000 0 0 1000);
ARRAY V[2,2] (1 0 0 1);
PARMS beta {5 2};
PARMS sigma_b {1 0 0 1};
PARMS var_e 1;
PRIOR beta ~ MVN(beta0, sigma0);
PRIOR Sigma_b ~ IWISH(2, V);
PRIOR var_e ~ IGAMMA(0.001, SCALE=0.001);
RANDOM b ~ MVN(beta, Sigma_b) SUBJECT=id;
mu = L + S * time;
MODEL y ~ NORMAL(mu, var=var_e);
RUN;

proc sgplot data=histnorm;
series x=Iteration y=beta1;
yaxis label="beta0";
run;

data _NULL_;
set histnorm;
file "C:\Documents and Settings\zzhang4\My Documents\My Dropbox\research\Bayesian flexible error\histnorm.txt";
put Iteration beta1 beta2 Sigma_b1 Sigma_b2 Sigma_b4 var_e;
run;


PROC MCMC DATA=ndata NMC=40000  NBI=5000 THIN=2
DIC SEED=17 INIT=RANDOM DIAG=(ESS GEWEKE(F1=.2 F2=.5)) 
STATISTICS(ALPHA=0.05)=(SUMMARY INTERVAL)
OUTPOST=histt; 
ODS SELECT PARAMETERS POSTSUMMARIES 
   TADPANEL GEWEKE ESS DIC POSTINTERVALS;
ARRAY b[2] L S;
ARRAY beta[2];
ARRAY Sigma_b[2,2];
ARRAY beta0[2] (0 0);
ARRAY sigma0[2,2] (1000 0 0 1000);
ARRAY V[2,2] (1 0 0 1);
PARMS beta {5 2};
PARMS sigma_b {1 0 0 1};
PARMS var_e 1;
PARMS df 5;
PRIOR beta ~ MVN(beta0, sigma0);
PRIOR Sigma_b ~ IWISH(2, V);
PRIOR var_e ~ IGAMMA(0.01, SCALE=0.01);
PRIOR df~UNIFORM(1,500);
RANDOM b ~ MVN(beta, Sigma_b) SUBJECT=id;
mu = L + S * time;
MODEL y ~ t(mu, var=var_e, df);
RUN;

proc sgplot data=histt;
series x=Iteration y=beta1;
yaxis label="beta0";
run;

data _NULL_;
set histt;
file "C:\Documents and Settings\zzhang4\My Documents\My Dropbox\research\Bayesian flexible error\histt.txt";
put Iteration beta1 beta2 Sigma_b1 Sigma_b2 Sigma_b4 var_e df;
run;