proc mcmc data=MFPT nbi=5000 nmc=50000 thin=10 monitor=(G littleg mu sigma); |
begincnst; J=4; c=1; endcnst; |
parms pi11 pi21 pi23; |
parms pi31 pi33 pi35 pi37; |
parms pi41 pi43 pi45 pi47 pi49 pi411 pi413 pi415; |
parms mu sigma / t(3); |
array p[]; |
p1 = pi11*pi21*pi31*pi41; p2 = pi11*pi21*pi31*(1-pi41); |
p3 = pi11*pi21*(1-pi31)*pi43; p4 = pi11*pi21*(1-pi31)*(1-pi43); |
p5 = pi11*(1-pi21)*pi33*pi45; p6 = pi11*(1-pi21)*pi33*(1-pi45); |
p7 = pi11*(1-pi21)*(1-pi33)*pi47; p8 = pi11*(1-pi21)*(1-pi33)*(1-pi47); |
p9 = (1-pi11)*pi23*pi35*pi49; p10 = (1-pi11)*pi23*pi35*(1-pi49); |
p11 = (1-pi11)*pi23*(1-pi35)*pi411; p12 = (1-pi11)*pi23*(1-pi35)*(1-pi411); |
p13 = (1-pi11)*(1-pi23)*pi37*pi413; p14 = (1-pi11)*(1-pi23)*pi37*(1-pi413); |
p15 = (1-pi11)*(1-pi23)*(1-pi37)*pi415; p16 = (1-pi11)*(1-pi23)*(1-pi37)*(1-pi415); |
prior pi11 ~ beta(c,c); prior pi21 pi23 ~ beta(c*2**2,c*2**2); |
prior pi31 pi33 pi35 pi37 ~ beta(c*3**2,c*3**2); |
prior pi41 pi43 pi45 pi47 pi49 pi411 pi413 pi415 ~ beta(c*4**2,c*4**2); |
prior mu ~ normal(50,sd=10); prior sigma ~ uniform(0,20); |
k = int(2**J * cdf("normal", y, mu, sigma) + 1); |
llike=J*log(2) + log(p[k]) + lpdfnorm(y, mu, sigma); |
model general(llike); |
array setID52; array G52; array littleg52; |
do i = 1 to dim(G) by 1; |
setID[i] = int(2**J * cdf("normal", 39 + i/2, mu, sigma) + 1); |
G[i] = p1*(setID[i] ge 2) + p2*(setID[i] ge 3) + p3*(setID[i] ge 4) |
+ p4*(setID[i] ge 5) + p5*(setID[i] ge 6) + p6*(setID[i] ge 7) |
+ p7*(setID[i] ge 8) + p8*(setID[i] ge 9) + p9*(setID[i] ge 10) |
+ p10*(setID[i] ge 11) + p11*(setID[i] ge 12) + p12*(setID[i] ge 13) |
+ p13*(setID[i] ge 14) + p14*(setID[i] ge 15) + p15*(setID[i] ge 16) |
+ p[setID{i}] *(2**J * cdf("normal", 39 + i/2, mu, sigma)- setID[i] + 1); |
littleg[i] = 2**J * p[setID{i}] * pdf("normal", 39 + i/2, mu, sigma); |
end; run; |