Documentation Example 19 for PROC MCMC
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: MCMCEX19 */
/* TITLE: Documentation Example 19 for PROC MCMC */
/* Implement a New Sampling Algorithm */
/* PRODUCT: STAT */
/* SYSTEM: ALL */
/* KEYS: */
/* PROCS: MCMC */
/* DATA: */
/* */
/* SUPPORT: Fang Chen */
/* REF: PROC MCMC, EXAMPLE 19 */
/* MISC: */
/****************************************************************/
title 'Implement a New Sampling Algorithm';
data inputdata;
input remiss cell smear infil li blast temp;
ind = _n_;
cnst = 1;
label remiss='Complete Remission';
datalines;
1 0.8 0.83 0.66 1.9 1.1 0.996
1 0.9 0.36 0.32 1.4 0.74 0.992
0 0.8 0.88 0.7 0.8 0.176 0.982
0 1 0.87 0.87 0.7 1.053 0.986
1 0.9 0.75 0.68 1.3 0.519 0.98
0 1 0.65 0.65 0.6 0.519 0.982
1 0.95 0.97 0.92 1 1.23 0.992
0 0.95 0.87 0.83 1.9 1.354 1.02
0 1 0.45 0.45 0.8 0.322 0.999
0 0.95 0.36 0.34 0.5 0 1.038
0 0.85 0.39 0.33 0.7 0.279 0.988
0 0.7 0.76 0.53 1.2 0.146 0.982
0 0.8 0.46 0.37 0.4 0.38 1.006
0 0.2 0.39 0.08 0.8 0.114 0.99
0 1 0.9 0.9 1.1 1.037 0.99
1 1 0.84 0.84 1.9 2.064 1.02
0 0.65 0.42 0.27 0.5 0.114 1.014
0 1 0.75 0.75 1 1.322 1.004
0 0.5 0.44 0.22 0.6 0.114 0.99
1 1 0.63 0.63 1.1 1.072 0.986
0 1 0.33 0.33 0.4 0.176 1.01
0 0.9 0.93 0.84 0.6 1.591 1.02
1 1 0.58 0.58 1 0.531 1.002
0 0.95 0.32 0.3 1.6 0.886 0.988
1 1 0.6 0.6 1.7 0.964 0.99
1 1 0.69 0.69 0.9 0.398 0.986
0 1 0.73 0.73 0.7 0.398 0.986
;
proc mcmc data=inputdata nmc=100000 propcov=quanew seed=17
outpost=mcmcout;
ods select PostSumInt ess;
parms beta0-beta6;
prior beta: ~ normal(0,var=25);
mu = beta0 + beta1*cell + beta2*smear +
beta3*infil + beta4*li + beta5*blast + beta6*temp;
p = cdf('normal', mu, 0, 1);
model remiss ~ bern(p);
run;
proc fcmp outlib=sasuser.funcs.uds;
/******************************************/
/* Generate left-truncated normal variate */
/******************************************/
function rltnorm(mu,sig,lwr);
if lwr<mu then do;
ans = lwr-1;
do while(ans<lwr);
ans = rand('normal',mu,sig);
end;
end;
else do;
mul = (lwr-mu)/sig;
alpha = (mul + sqrt(mul**2 + 4))/2;
accept=0;
do while(accept=0);
z = mul + rand('exponential')/alpha;
lrho = -(z-alpha)**2/2;
u = rand('uniform');
lu = log(u);
if lu <= lrho then accept=1;
end;
ans = sig*z + mu;
end;
return(ans);
endsub;
/*******************************************/
/* Generate right-truncated normal variate */
/*******************************************/
function rrtnorm(mu,sig,uppr);
ans = 2*mu - rltnorm(mu,sig, 2*mu-uppr);
return(ans);
endsub;
run;
/* define the HH algorithm in PROC FCMP. */
%let n = 27;
%let p = 7;
options cmplib=sasuser.funcs;
proc fcmp outlib=sasuser.funcs.uds;
subroutine HH(beta[*],Z[*],B[*],x[*,*],y[*],W[*],sQ[*],S[*,*],L[*,*]);
outargs beta, Z, B;
array T[&p] / nosym;
do j=1 to &n;
zold = Z[j];
m = 0;
do k = 1 to &p;
m = m + X[j,k] * B[k];
end;
m = m - W[j]*(Z[j]-m);
if (y[j]=1) then
Z[j] = rltnorm(m,sQ[j],0);
else
Z[j] = rrtnorm(m,sQ[j],0);
diff = Z[j] - zold;
do k = 1 to &p;
B[k] = B[k] + diff * S[k,j];
end;
end;
do j=1 to &p;
T[j] = rand('normal');
end;
call mult(L,T,T);
call addmatrix(B,T,beta);
endsub;
run;
options cmplib=sasuser.funcs;
proc mcmc data=inputdata nmc=5000 monitor=(beta) outpost=hhout;
ods select PostSumInt ess;
array xtx[&p,&p]; /* work space */
array xt[&p,&n]; /* work space */
array v[&p,&p]; /* work space */
array HatMat[&n,&n]; /* work space */
array S[&p,&n]; /* V * Xt */
array W[&n];
array y[1]/ nosymbols; /* y stores the response variable */
array x[1]/ nosymbols; /* x stores the explanatory variables */
array sQ[&n]; /* sqrt of the diagonal elements of Q */
array B[&p]; /* conditional mean of beta */
array L[&p,&p]; /* Cholesky decomp of conditional cov */
array Z[&n]; /* latent variables Z */
array beta[&p] beta0-beta6; /* regression coefficients */
begincnst;
call streaminit(83101);
if ind=1 then do;
rc = read_array("inputdata", x, "cnst", "cell", "smear", "infil",
"li", "blast", "temp");
rc = read_array("inputdata", y, "remiss");
call identity(v);
call mult(v, 25, v);
call transpose(x,xt);
call mult(xt,x,xtx);
call inv(v,v);
call addmatrix(xtx,v,xtx);
call inv(xtx,v);
call chol(v,L);
call mult(v,xt,S);
call mult(x,S,HatMat);
do j=1 to &n;
H = HatMat[j,j];
W[j] = H/(1-H);
sQ[j] = sqrt(W[j] + 1);
end;
do j=1 to &n;
if(y[j]=1) then
Z[j] = rltnorm(0,1,0);
else
Z[j] = rrtnorm(0,1,0);
end;
call mult(S,Z,B);
end;
endcnst;
uds HH(beta,Z,B,x,y,W,sQ,S,L);
parms z: beta: 0 B1-B7 / uds;
prior z: beta: B1-B7 ~ general(0);
model general(0);
run;