### R code from vignette source 'mcmc_practice.Rnw' ################################################### ### code chunk number 1: pima ################################################### library(MASS) data(Pima.te) m0=glm(type~-1+glu+bp,family=binomial(link="probit"),data=Pima.te,x=TRUE) m1=glm(type~-1+glu+bp+ped,family=binomial(link="probit"),data=Pima.te,x=TRUE) ################################################### ### code chunk number 2: mcmc_practice.Rnw:168-173 ################################################### y=m0$y lik0=function(beta) exp(sum ( y*pnorm(m0$x%*%beta,log.p=TRUE)+ (1-y)*pnorm(m0$x%*%beta,log.p=TRUE,lower.tail=FALSE)) ) lik1=function(beta) exp(sum ( y*pnorm(m1$x%*%beta,log.p=TRUE)+ (1-y)*pnorm(m1$x%*%beta,lower.tail=FALSE,log.p=TRUE) )) ################################################### ### code chunk number 3: mcmc_practice.Rnw:180-189 ################################################### library(mvtnorm) S0=solve(t(m0$x)%*%m0$x)*length(y) S1=solve(t(m1$x)%*%m1$x)*length(y) prior0.sim=rmvnorm(10000,sigma=S0) prior1.sim=rmvnorm(10000,sigma=S1) marg0=apply(prior0.sim,FUN=lik0,MAR=1) marg1=apply(prior1.sim,FUN=lik1,MAR=1) B01=sum(marg0)/sum(marg1) B01 ################################################### ### code chunk number 4: mcmc_practice.Rnw:204-218 ################################################### S0.g=summary(m0)$cov.unscaled S1.g=summary(m1)$cov.unscaled g0.sim=rmvnorm(10000,mean=coef(m0), sigma=S0.g) g1.sim=rmvnorm(10000,mean=coef(m1), sigma=S1.g) importance0=function(beta) dmvnorm(beta,mean=c(0,0),sigma=S0)/ dmvnorm(beta,mean=coef(m0),sigma=S0.g) importance1=function(beta) dmvnorm(beta,mean=c(0,0,0),sigma=S1)/ dmvnorm(beta,mean=coef(m1),sigma=S1.g) f0=function(beta) lik0(beta)*importance0(beta) f1=function(beta) lik1(beta)*importance1(beta) marg0=apply(g0.sim,FUN=f0,MAR=1) marg1=apply(g1.sim,FUN=f1,MAR=1) B01=sum(marg0)/sum(marg1) B01 ################################################### ### code chunk number 5: mcmc_practice.Rnw:331-347 ################################################### x=read.table("caterpillar.dat") y=log(x[,11]) x=x[,-11] n=length(y) one=rep(1,n); x=cbind(one,x) x=as.matrix(x) cost=100 p=ncol(x) ss=sum(y^2) marg.gamma=function(gam){ gam=as.logical(gam) xg=x[,gam] marg=((cost+1)^(-(sum(gam)+1)/2))*(ss-cost/(cost+1)* t(y)%*%xg%*%solve(t(xg)%*%xg) %*%t(xg) %*%y)^(-n/2) marg } ################################################### ### code chunk number 6: mcmc_practice.Rnw:355-371 ################################################### gibbs=function(nsim,gamma.init){ gam=gamma.init out=matrix(nrow=nsim,ncol=length(gam)) for (iter in 1:nsim){ for (j in 2:p){ gam[j]=1 p1=marg.gamma(gam) gam[j]=0 p0=marg.gamma(gam) gam[j]=sample(x=c(0,1),prob=c(p0,p1),size=1)} out[iter,]=gam} out} gamma.init=c(1,sample(c(0,1),size=(ncol(x)-1),rep=TRUE)) out=gibbs(nsim=20000,gamma.init) sort(table(apply(out,FUN=paste,MAR=1,collapse="")), decreasing=TRUE)[1:10]/nrow(out) ################################################### ### code chunk number 7: mh1 ################################################### temp=c(66,70,69,68,67,72,73,70,57,63,70,78,67,53,67,75,70,81,76,79,75,76,58) failure=c(0,1,0,0,0,0,0,0,1,1,1,0,0,2,0,0,0,0,0,0,2,0,1) loglik <- function(beta) {sum(failure * (beta[1] + beta[2] * temp)) -sum(6 * log(1 + exp(beta[1] + beta[2] *temp)))} beta0=seq(-5,15,length=200) beta1=seq(-0.25,0.05,length=200) l.val=apply(expand.grid(beta0,beta1),FUN=loglik,MAR=1) l.val=matrix(l.val,ncol=length(beta1),nrow=length(beta0)) contour(beta0,beta1,exp(l.val-max(l.val)),levels=c(0.01,0.1,0.5,0.9)) ################################################### ### code chunk number 8: mcmc_practice.Rnw:485-495 ################################################### mh <- function(nsim, s1, s2, b.init) { mh.out <- matrix(ncol = 2, nrow = nsim) b <- b.init for (i in 1:nsim) { b.p <- c(rnorm(1, b[1], s1), rnorm(1,b[2], s2)) if (runif(1) < exp(loglik(b.p) - loglik(b))) b <- b.p mh.out[i, ] <- b } mh.out} ################################################### ### code chunk number 9: mh2 ################################################### mh.out <- mh(nsim = 5000, s1 = 1, s2 = 1, b.init = c(0,0)) par(mfrow = c(2, 2)) plot(mh.out[, 1],type="l") plot(mh.out[, 2],type="l") contour(beta0,beta1,exp(l.val-max(l.val)),levels=c(0.01,0.1,0.5,0.9)) points(mh.out,col=2) ################################################### ### code chunk number 10: mh3 ################################################### mh.out=mh(nsim = 5000, s1 = 2, s2 = 0.1, b.init = c(0,0)) par(mfrow = c(2, 2)) plot(mh.out[, 1],type="l") plot(mh.out[, 2],type="l") contour(beta0,beta1,exp(l.val-max(l.val)),levels=c(0.01,0.1,0.5,0.9)) points(mh.out,col=2) ################################################### ### code chunk number 11: mh4 ################################################### par(mfrow = c(2, 2)) mh.out=mh(nsim = 10000, s1 = sd(mh.out[,1]), s2 = sd(mh.out[,2]), b.init = c(0,0)) plot(mh.out[, 1],type="l") plot(mh.out[, 2],type="l") contour(beta0,beta1,exp(l.val-max(l.val)),levels=c(0.01,0.1,0.5,0.9)) points(mh.out,col=2) ################################################### ### code chunk number 12: mcmc_practice.Rnw:591-602 ################################################### library(mvtnorm) mhd <- function(nsim, V, b.init) { mh.out <- matrix(ncol = 2, nrow = nsim) b <- b.init for (i in 1:nsim) { b.p <- rmvnorm(n = 1, mean = b, sigma = V) if (runif(1) < exp(loglik(b.p) - loglik(b))) b <- b.p mh.out[i, ] <- b } mh.out} ################################################### ### code chunk number 13: mh5 ################################################### par(mfrow=c(2,2)) v <- var(mh.out) mh.out <- mhd(nsim = 10000, V = v, b.init = c(0,0)) plot(mh.out[, 1],type="l") plot(mh.out[, 2],type="l") contour(beta0,beta1,exp(l.val-max(l.val)),levels=c(0.01,0.1,0.5,0.9)) points(mh.out,col=2) ################################################### ### code chunk number 14: mh6 ################################################### apply(mh.out,FUN=mean,MAR=2) apply(mh.out,FUN=sd,MAR=2) par(mfrow = c(1, 2)) plot(density(mh.out[, 1], bw = 0.5), xlab = "beta0",main="") plot(density(mh.out[, 2], bw = 0.01), xlab = "beta1",main="") ################################################### ### code chunk number 15: predmh ################################################### pred=rbinom(nrow(mh.out),size=6, prob=exp(mh.out%*%(c(1,31)))/(1+exp(mh.out%*%(c(1,31))))) mean(pred) barplot(table(pred)/nrow(mh.out),col="sienna2") ################################################### ### code chunk number 16: mcmc_practice.Rnw:707-713 ################################################### b <- function(p) { b1 <- (1/20) * log(p[2] * (1 - p[1])/(p[1] * (1 - p[2]))) b2 <- log(p[2]/(1 - p[2])) - 60 * b1 return(c(b2, b1)) } ################################################### ### code chunk number 17: mcmc_practice.Rnw:763-771 ################################################### p.sim <- matrix(rbeta(2000, 1, 1), ncol = 2) b.sim <- t(apply(p.sim, MAR = 1, FUN = b)) sum(p.sim[,1]*exp(apply(b.sim, FUN = loglik, MAR = 1)))/ sum(exp(apply(b.sim, FUN = loglik, MAR = 1))) p.sim <- matrix(rbeta(2000, 1, 1), ncol = 2) b.sim <- t(apply(p.sim, MAR = 1, FUN = b)) sum(p.sim[,1]*exp(apply(b.sim, FUN = loglik, MAR = 1)))/ sum(exp(apply(b.sim, FUN = loglik, MAR = 1))) ################################################### ### code chunk number 18: mcmc_practice.Rnw:784-794 ################################################### est=c() for (i in 1:1000){ p.sim <- matrix(rbeta(2000, 1, 1), ncol = 2) b.sim <- t(apply(p.sim, MAR = 1, FUN = b)) est[i]=sum(p.sim[,1]*exp(apply(b.sim, FUN = loglik, MAR = 1)))/ sum(exp(apply(b.sim, FUN = loglik, MAR = 1)))} mean(est) sd(est) ################################################### ### code chunk number 19: binomial2 ################################################### plot(b.sim,xlab=expression(beta[0]),ylab=expression(beta[1]), cex=0.75) contour(beta0,beta1,exp(l.val-max(l.val)), levels=c(0.01,0.1,0.5,0.9), add=TRUE,col=2,lwd=2) ################################################### ### code chunk number 20: binomial3 ################################################### w <- exp(apply(b.sim, FUN = loglik, MAR = 1)) s <- sample(x = 1000, size = 1000, prob = w, rep = T) plot(b.sim[s,],xlab=expression(beta[0]),ylab=expression(beta[1])) contour(beta0,beta1,exp(l.val-max(l.val)),levels=c(0.01,0.1,0.5,0.9), add=TRUE,col=2,lwd=2) ################################################### ### code chunk number 21: mcmc_practice.Rnw:894-907 ################################################### m1=c(mean(p.sim[s, 1]), mean(p.sim[s, 1]^2)) m2=c(mean(p.sim[s, 2]), mean(p.sim[s, 2]^2)) g1.sim=rbeta(1000,shape1=m1[1]*(m1[2]-m1[1])/(m1[1]^2-m1[2]), shape2=(1-m1[1])*(m1[2]-m1[1])/(m1[1]^2-m1[2])) g2.sim=rbeta(1000,shape1=m2[1]*(m2[2]-m2[1])/(m2[1]^2-m2[2]), shape2=(1-m2[1])*(m2[2]-m2[1])/(m2[1]^2-m2[2])) g.sim=cbind(g1.sim,g2.sim); b.sim=t(apply(g.sim,FUN=b,MAR=1)) prop1.g=dbeta(g1.sim,shape1=m1[1]*(m1[2]-m1[1])/(m1[1]^2-m1[2]), shape2=(1-m1[1])*(m1[2]-m1[1])/(m1[1]^2-m1[2])) prop2.g=dbeta(g2.sim,shape1=m2[1]*(m2[2]-m2[1])/(m2[1]^2-m2[2]), shape2=(1-m2[1])*(m2[2]-m2[1])/(m2[1]^2-m2[2])) w=exp(apply(b.sim,FUN=loglik,MAR=1))/(prop1.g*prop2.g) sum(g1.sim*w)/sum(w) ################################################### ### code chunk number 22: mcmc_practice.Rnw:918-934 ################################################### est=c() for (i in 1:1000){ g1.sim=rbeta(1000,shape1=m1[1]*(m1[2]-m1[1])/(m1[1]^2-m1[2]), shape2=(1-m1[1])*(m1[2]-m1[1])/(m1[1]^2-m1[2])) g2.sim=rbeta(1000,shape1=m2[1]*(m2[2]-m2[1])/(m2[1]^2-m2[2]), shape2=(1-m2[1])*(m2[2]-m2[1])/(m2[1]^2-m2[2])) g.sim=cbind(g1.sim,g2.sim); b.sim=t(apply(g.sim,FUN=b,MAR=1)) prop1.g=dbeta(g1.sim,shape1=m1[1]*(m1[2]-m1[1])/(m1[1]^2-m1[2]), shape2=(1-m1[1])*(m1[2]-m1[1])/(m1[1]^2-m1[2])) prop2.g=dbeta(g2.sim,shape1=m2[1]*(m2[2]-m2[1])/(m2[1]^2-m2[2]), shape2=(1-m2[1])*(m2[2]-m2[1])/(m2[1]^2-m2[2])) w=exp(apply(b.sim,FUN=loglik,MAR=1))/(prop1.g*prop2.g) est[i]=sum(g1.sim*w)/sum(w)} mean(est) sd(est) ################################################### ### code chunk number 23: binomial4 ################################################### plot(b.sim,xlab=expression(beta[0]),ylab=expression(beta[1])) contour(beta0,beta1,exp(l.val-max(l.val)), levels=c(0.01,0.1,0.5,0.9), add=TRUE,col=2,lwd=2) ################################################### ### code chunk number 24: binomial5 ################################################### s <- sample(x = 1000, size = 1000, prob = w, rep = T) par(mfrow=c(1,3)) plot(b.sim[s,],xlab=expression(beta[0]),ylab=expression(beta[1])) contour(beta0,beta1,exp(l.val-max(l.val)),levels=c(0.01,0.1,0.5,0.9), add=TRUE,col=2,lwd=2) plot(density(b.sim[s,1]),xlab=expression(beta[0]),cex.lab=1.25, main="") plot(density(b.sim[s,2]),xlab=expression(beta[1]),cex.lab=1.25, main="") ################################################### ### code chunk number 25: mcmc_practice.Rnw:1160-1162 ################################################### library(boot) data(coal) ################################################### ### code chunk number 26: coal1 ################################################### anni = seq(1851, 1963, 1) y = hist(coal$date, breaks = anni)$counts ################################################### ### code chunk number 27: mcmc_practice.Rnw:1195-1218 ################################################### kpost.cond = function(y, mean1, mean2) { exp((mean1 - mean2) * seq(1:length(y))) * exp(log(mean2/mean1) * cumsum(y))} cppost = function(y, nsim, a1, a2, c1, c2, d1, d2) { sim = matrix(nrow = nsim, ncol = 5) lambda = mean(y) theta = mean(y) for (j in 1:nsim) { k = sample(x = length(y), size = 1, prob = kpost.cond(y, mean1 = lambda, mean2 = theta)) b1 = rgamma(n = 1, shape = a1 + c1, rate = theta + d1) b2 = rgamma(n = 1, shape = a2 + c2, rate = lambda + d2) theta = rgamma(n = 1, shape = a1 + cumsum(y)[k], rate = b1 + k) lambda = rgamma(n = 1, shape = a2 + (sum(y) - cumsum(y)[k]), rate = b2 + (length(y) - k)) sim[j, ] = c(k, theta, lambda, b1, b2) } return(sim)} ################################################### ### code chunk number 28: coal2 ################################################### out = cppost(y = y, nsim = 2000, a1 = 1, a2 = 1,c1 = 1, c2 = 1, d1 = 1, d2 = 1) par(mfrow = c(2, 3)) par(mar = c(3, 3, 3, 3)) plot(out[, 1], ylab = "k", col = 2, cex = 1.2) plot(out[, 2], ylab = expression(theta), type = "l",col = 2, cex = 1.2) plot(out[, 3], ylab = expression(lambda), type = "l",col = 2, cex = 1.2) barplot(table(out[100:2000, 1] + 1851), col = "orange") plot(density(out[100:2000, 2]), main = "", xlab = expression(theta),cex = 1.2) plot(density(out[100:2000, 3]), main = "", xlab = expression(lambda),cex = 1.2) ################################################### ### code chunk number 29: samplemix ################################################### set.seed(7) mu1 = 2.5; mu2 = 0; p = 0.7; n = 500; pbar = 1 - p u = runif(n) sampl = rnorm(n) + (u <= p) * mu1 + (u > p) * mu2 hist(sampl,nclass=20,prob=TRUE,main="",col="orange") lines(density(sampl),col="purple",lwd=2) ################################################### ### code chunk number 30: mix2 ################################################### mu1 = mu2 = seq(min(sampl), max(sampl), 0.1) mo1 = mu1 %*% t(rep(1, length(mu2))) mo2 = rep(1, length(mu2)) %*% t(mu2) ca1 = -0.5 * mo1 * mo1;ca2 = -0.5 * mo2 * mo2 like = 0 * mo1 for (i in 1:n) like = like + log(p * exp(ca1 + sampl[i] * mo1) + pbar * exp(ca2 + sampl[i] * mo2)) image(mu1,mu2,like,col=heat.colors(250)) contour(mu1,mu2,like,add=TRUE,levels=seq(min(like),max(like),length=50)) ################################################### ### code chunk number 31: mcmc_practice.Rnw:1482-1504 ################################################### gibbsmean=function (p, datha, niter = 10^4) { n = length(datha) z = rep(0, n) ssiz = nxj = rep(0, 2) mug = matrix(0, nrow = niter + 1, ncol = 2) mug[1, ] = rep(mean(datha), 2) for (i in 2:(niter + 1)) { for (t in 1:n) { prob = c(p, 1 - p) * dnorm(datha[t], mean = mug[i - 1, ]) z[t] = sample(c(1, 2), size = 1, prob = prob) } for (j in 1:2) { ssiz[j] = 1 + sum(z == j) nxj[j] = sum(as.numeric(z == j) * datha) } mug[i, ] = rnorm(2, mean = nxj/ssiz, sd = sqrt(1/ssiz)) } mug } ################################################### ### code chunk number 32: mix3 ################################################### simu1=gibbsmean(0.7,sampl,niter=10^3) simu2=gibbsmean(0.7,sampl,niter=10^3) image(mu1,mu2,like,col=heat.colors(250)) contour(mu1,mu2,like,add=TRUE,levels=seq(min(like),max(like),length=50)) points(simu1,col=2); points(simu2,col=3) ################################################### ### code chunk number 33: license1 ################################################### library(bayess) data(datha) hist(as.matrix(datha),nclass=200,prob=TRUE,main="") ################################################### ### code chunk number 34: mcmc_practice.Rnw:1690-1717 ################################################### gibbsknorm=function (niter, dat, mix) { n = length(dat);k = mix$k z = rep(0, n);ssiz = nxj = ssum = rep(0, k) mug = sigg = prog = matrix(0, nrow = niter, ncol = k) prog[1, ] = rep(1, k)/k; mug[1, ] = rep(mix$mu, k) sigg[1, ] = rep(mix$sig, k) for (i in 1:(niter - 1)) { for (t in 1:n) { prob = prog[i, ] * dnorm(dat[t], mean = mug[i, ], sd = sqrt(sigg[i, ])) if (sum(prob) == 0) prob = rep(1, k)/k z[t] = sample(1:k, 1, prob = prob)} for (j in 1:k) { ssiz[j] = sum(z == j) nxj[j] = sum(as.numeric(z == j) * dat)} mug[i + 1, ] = rnorm(k, (mean(dat) + nxj)/(ssiz + 1), sqrt(sigg[i, ]/(ssiz + 1))) for (j in 1:k) ssum[j] = sum(as.numeric(z == j) * (dat - nxj[j]/ssiz[j])^2) mxj=nxj/ssiz; mxj[ssiz==0]=0; ssum[ssiz==0]=0 sigg[i + 1, ] = 1/rgamma(k, shape = 0.5 * (20 + ssiz), rate = var(dat) + (0.5 * ssum + 0.5 * ssiz/(ssiz + 1) * (mean(dat) - mxj)^2)) prog[i + 1, ] = rdirichlet(1, par = ssiz + 0.5) } list(k = k, mu = mug, sig = sigg, p = prog) } ################################################### ### code chunk number 35: license2 ################################################### datha=as.matrix(datha) mix=list(k=3,mu=mean(datha),sig=var(datha)) simu=gibbsknorm(1000,dat=as.matrix(datha),mix=mix) hist(datha,prob=TRUE,xlab="",ylab="",nclass=100,main="") x=y=seq(min(datha),max(datha),length=150) yy=matrix(0,ncol=150,nrow=1000) for (i in 1:150){ yy[,i]=rowSums(simu$p*dnorm(x[i], mean=simu$mu,sd=sqrt(simu$sig))) y[i]=mean(yy[,i])} for (t in 501:1000) lines(x,yy[t,],col="gold") lines(x,y,lwd=2.5,col="sienna2") ################################################### ### code chunk number 36: license3 ################################################### par(mfrow=c(3,2)) par(mar=c(2,2,2,2)) plot(simu$mu[,2],ylim=range(simu$mu), ylab=expression(mu[i]),xlab="n",type="l",col="sienna3") lines(simu$mu[,1],col="gold4") lines(simu$mu[,3],col="steelblue") plot(density(simu$mu[-(1:10),2]),xlim=range(simu$mu), xlab=expression(mu[i]),type="l",col="sienna3",ylab="",main="",lwd=3) lines(density(simu$mu[-(1:10),1]),col="gold4",lwd=3) lines(density(simu$mu[-(1:10),3]),col="steelblue",lwd=3) plot(simu$p[,2],ylim=range(simu$p), ylab=expression(p[i]),xlab="n",type="l",col="sienna3") lines(simu$p[,1],col="gold4") lines(simu$p[,3],col="steelblue") plot(density(simu$p[-(1:10),2]),xlim=range(simu$p), xlab=expression(p[i]),ylab="",type="l",col="sienna3",main="",lwd=3) lines(density(simu$p[-(1:10),1]),col="gold4",lwd=3) lines(density(simu$p[-(1:10),3]),col="steelblue",lwd=3) plot(simu$sig[,2],ylim=range(simu$sig), ylab=expression(sigma^2[i]),xlab="n",type="l",col="sienna3") lines(simu$sig[,1],col="gold4") lines(simu$sig[,3],col="steelblue") plot(density(simu$sig[-(1:10),2]),xlim=range(simu$sig), xlab=expression(sig[i]),ylab="",type="l",col="sienna3",main="",lwd=3) lines(density(simu$sig[-(1:10),1]),col="gold4",lwd=3) lines(density(simu$sig[-(1:10),3]),col="steelblue",lwd=3) ################################################### ### code chunk number 37: mcmc_practice.Rnw:1831-1833 ################################################### dat.small=as.matrix(datha)[1:150,] simu=gibbsknorm(10000,dat=dat.small,mix=mix) ################################################### ### code chunk number 38: license4 ################################################### par(mfrow=c(3,2)) par(mar=c(2,2,2,2)) plot(simu$mu[,2],ylim=range(simu$mu), ylab=expression(mu[i]),xlab="n",type="l",col="sienna3") lines(simu$mu[,1],col="gold4") lines(simu$mu[,3],col="steelblue") plot(density(simu$mu[-(1:10),2]),xlim=range(simu$mu), xlab=expression(mu[i]),type="l",col="sienna3",ylab="",main="",lwd=3) lines(density(simu$mu[-(1:10),1]),col="gold4",lwd=3) lines(density(simu$mu[-(1:10),3]),col="steelblue",lwd=3) plot(simu$p[,2],ylim=range(simu$p), ylab=expression(p[i]),xlab="n",type="l",col="sienna3") lines(simu$p[,1],col="gold4") lines(simu$p[,3],col="steelblue") plot(density(simu$p[-(1:10),2]),xlim=range(simu$p), xlab=expression(p[i]),ylab="",type="l",col="sienna3",main="",lwd=3) lines(density(simu$p[-(1:10),1]),col="gold4",lwd=3) lines(density(simu$p[-(1:10),3]),col="steelblue",lwd=3) plot(simu$sig[,2],ylim=range(simu$sig), ylab=expression(sigma^2[i]),xlab="n",type="l",col="sienna3") lines(simu$sig[,1],col="gold4") lines(simu$sig[,3],col="steelblue") plot(density(simu$sig[-(1:10),2]),xlim=range(simu$sig), xlab=expression(sig[i]),ylab="",type="l",col="sienna3",main="",lwd=3) lines(density(simu$sig[-(1:10),1]),col="gold4",lwd=3) lines(density(simu$sig[-(1:10),3]),col="steelblue",lwd=3) ################################################### ### code chunk number 39: mcmc_practice.Rnw:1916-1928 ################################################### log.post=function(simu,data){ n=length(data) mu=simu$mu;sig=simu$sig;p=simu$p;k=simu$k logpost=rep(0,nrow(mu)) for (it in 1:nrow(mu)){ lik=matrix(0,n,k) for(j in 1:k) lik[,j]=p[it,j]*dnorm(x=data,mean=mu[it,j],sd=sqrt(sig[it,j])) logpost[it]=sum(log(rowSums(lik)))+sum(dnorm(mu[it,],mean(data),sqrt(sig[it,]),log=TRUE))- (10+1)*sum(log(sig[it,]))-sum(var(data)/sig[it,])+0.5*sum(log(p[it,])) } logpost} lp=log.post(simu,dat.small) ################################################### ### code chunk number 40: mcmc_practice.Rnw:1939-1941 ################################################### ind.map=order(lp,decreasing=TRUE)[1] map=list(mu=simu$mu[ind.map,],sig=simu$sig[ind.map,],p=simu$p[ind.map,]) ################################################### ### code chunk number 41: mcmc_practice.Rnw:1946-1950 ################################################### pz.x.star=matrix(0,length(dat.small),3) for (i in 1:length(dat.small)){ pz.x.star[i,]=map$p*dnorm(dat.small[i],mean=map$mu,sd=sqrt(map$sig)) pz.x.star[i,]=pz.x.star[i,]/sum(pz.x.star[i,])} ################################################### ### code chunk number 42: pivoting ################################################### or.mu=or.sig=or.p=matrix(0,ncol=3,nrow=1000) pz.x=matrix(0,length(dat.small),3) library(combinat) perma=permn(3) for ( it in 1:1000){ kl.dist=rep(0,factorial(3)) for (i in 1:length(dat.small)){ pz.x[i,]=simu$p[it]*dnorm(dat.small[i],mean=simu$mu[it,],sd=sqrt(simu$sig[it,])) pz.x[i,]=pz.x[i,]/sum(pz.x[i,]) for (tau in 1:factorial(3)) kl.dist[tau]=kl.dist[tau]+sum(pz.x.star[i,]*log(pz.x[i,perma[[tau]]]))} best=order(kl.dist,decreasing=TRUE)[1] or.mu[it,]=simu$mu[it,perma[[best]]] or.sig[it,]=simu$sig[it,perma[[best]]] or.p[it,]=simu$p[it,perma[[best]]] } ################################################### ### code chunk number 43: mcmc_practice.Rnw:1991-1991 ################################################### ################################################### ### code chunk number 44: license5 ################################################### par(mfrow=c(3,2)) par(mar=c(2,2,2,2)) simu$mu=or.mu;simu$sig=or.sig;simu$p=or.p plot(simu$mu[,2],ylim=range(simu$mu), ylab=expression(mu[i]),xlab="n",type="l",col="sienna3") lines(simu$mu[,1],col="gold4") lines(simu$mu[,3],col="steelblue") plot(density(simu$mu[-(1:10),2]),xlim=range(simu$mu), xlab=expression(mu[i]),type="l",col="sienna3",ylab="",main="",lwd=3) lines(density(simu$mu[-(1:10),1]),col="gold4",lwd=3) lines(density(simu$mu[-(1:10),3]),col="steelblue",lwd=3) plot(simu$p[,2],ylim=range(simu$p), ylab=expression(p[i]),xlab="n",type="l",col="sienna3") lines(simu$p[,1],col="gold4") lines(simu$p[,3],col="steelblue") plot(density(simu$p[-(1:10),2]),xlim=range(simu$p), xlab=expression(p[i]),ylab="",type="l",col="sienna3",main="",lwd=3) lines(density(simu$p[-(1:10),1]),col="gold4",lwd=3) lines(density(simu$p[-(1:10),3]),col="steelblue",lwd=3) plot(simu$sig[,2],ylim=range(simu$sig), ylab=expression(sigma^2[i]),xlab="n",type="l",col="sienna3") lines(simu$sig[,1],col="gold4") lines(simu$sig[,3],col="steelblue") plot(density(simu$sig[-(1:10),2]),xlim=range(simu$sig), xlab=expression(sig[i]),ylab="",type="l",col="sienna3",main="",lwd=3) lines(density(simu$sig[-(1:10),1]),col="gold4",lwd=3) lines(density(simu$sig[-(1:10),3]),col="steelblue",lwd=3) ################################################### ### code chunk number 45: mcmc_practice.Rnw:2112-2115 ################################################### library(SMPracticals) data(sticky) head(sticky,7) ################################################### ### code chunk number 46: mcmc_practice.Rnw:2121-2126 ################################################### library(nlme) gsticky=groupedData(y~1|subject,data=sticky) m=lme(y~1,data=sticky,random= ~1|subject) m c(m$coef$random$subject+m$coef$fixed) ################################################### ### code chunk number 47: mcmc_practice.Rnw:2132-2151 ################################################### gibbs.re=function(nsim,hyper,data){ T=nlevels(data$subject); R=nrow(data)/T mu0=hyper[1];tau2=hyper[2] alpha=hyper[3]; beta=hyper[4] a.the=hyper[5]; b.the=hyper[6] ybar=tapply(data$y,FUN=mean,INDEX=data$subject) theta=ybar; mu=mean(theta) out=matrix(nrow=nsim,ncol=T+3) for ( iter in 1:nsim){ s2.the=1/rgamma(1,shape=a.the+T/2,rate=b.the+sum((theta-mu)^2)/2) s2=1/rgamma(1,shape=alpha+R*T/2,rate=beta+sum((data$y-rep(theta,R))^2)/2) mean.mu=(mean(theta)*T/s2.the+mu0/tau2)/(T/s2.the+1/tau2) var.mu=(T/s2.the+1/tau2)^(-1) mu=rnorm(1,mean.mu,sqrt(var.mu)) mean.the=(ybar*R/s2+mu/s2.the)/(R/s2+1/s2.the) var.the=(R/s2+1/s2.the)^(-1) theta=rnorm(T,mean=mean.the,sd=sqrt(var.the)) out[iter,]=c(s2.the,s2,mu,theta)} return(out) } ################################################### ### code chunk number 48: mcmc_practice.Rnw:2159-2166 ################################################### out=gibbs.re(nsim=10000,hyper=c(0,1000,0.5,1,0.5,1),data=sticky) out=out[1000:10000,] m=colMeans(out) sqrt(m[1:2]) m[3] m[4:9] apply(out,FUN=sd,MAR=2) ################################################### ### code chunk number 49: mcmc_practice.Rnw:2176-2177 ################################################### quantile(out[,1]/out[,2],c(0.05,0.95)) ################################################### ### code chunk number 50: mcmc_practice.Rnw:2234-2249 ################################################### y=rt(50,df=5) gibbs.t=function(nsim,y,nu.init){ out=double(nsim) n=length(y) nu=nu.init for (i in 1:nsim){ z=rgamma(n,shape=(nu+1)/2,rate=(nu+y^2)/2) nu.prop=rnorm(1,nu,sd=1) if (nu.prop<1 | nu.prop>100) nu.prop=nu num=sum(dgamma(z,nu.prop/2,nu.prop/2,log=TRUE)) den=sum(dgamma(z,nu/2,nu/2,log=TRUE)) if (runif(1)