Simple S-Plus code for Gibbs sampling analysis of a two-component normal
mixture with known variance and one component centred at zero, the other
having a positive mean.
# Simulate a simple normal mixture, one mean zero and known variance
# Plot resulting histogram and overlay density curve
- n_200
- mux_2.5
- sigma_1
- probx_0.75
- zx_runif(n)>probx
- y_sort(rnorm(n,zx*mux,sigma))
- hist(y,nclass=20,probability=T,col=3)
- x_seq(-3*sigma,mux+3*sigma,length=100)
- lines(x,probx*dnorm(x,0,sigma)+(1-probx)*dnorm(x,mux,sigma))
# Gibbs sampling for inference in above problem
- # MCMC sample size nmc; some initial "burn-in" iterations to discard, nit
- nmc_2000; nit_200
- # initial values
- mu_5; mu.s_NULL
- prob_.5; prob.s_NULL
- z_c(rep(0,n/2),rep(1,n/2))
- pz_matrix(0,n,nmc+nit)
- for (i in 1:(nmc+nit)) {
- # simulate mu and prob given z
-
- tot_sum(z)
-
- m_sum(z*y)/tot
-
- v_sigma/sqrt(tot)
- # simulate mu and pi given z (note the "tricky" bit for truncated normal
-
- a_pnorm(-m/v); b_pnorm((5-m)/v)
-
- mu_m+v*qnorm(a+runif(1)*(b-a))
-
- prob_rbeta(1,n-tot+1,tot+1)
- # then simulate new z given mu and prob (note vector form)
-
- prob.s_c(prob.s,prob); mu.s_c(mu.s,mu)
-
- qz_(1-prob)/(1-prob+prob*exp(-mu*(y-mu/2)/sigma^2))
-
- pz[,i]_qz
-
- z_runif(n)<=qz
-
- }
# plot MCMC iterations at start to check convergence
postscript(file="mix.ps",horizontal=T)
- par(mfrow=c(2,1),bty='n',pty='r',oma=c(1,1,1,1))
- tsplot(mu.s[1:(2*nit)],xlab="MCMC iteration",ylab="mu")
- tsplot(prob.s[1:(2*nit)],xlab="MCMC iteration",ylab="pi")
# Drop first nit draws for "burn-in" of MCMC and summarise the rest
- mu.s_mu.s[-(1:nit)]
- prob.s_prob.s[-(1:nit)]
- pz_pz[,-(1:nit)]%*%rep(1,nmc)/nmc
# histograms of posterior samples
- postscript(file="mix.ps",horizontal=T)
- par(mfrow=c(2,2),bty='n')
- hist(y,nclass=20,probability=T,xlab='data y')
- x_seq(-3*sigma,mux+3*sigma,length=100)
- lines(x,probx*dnorm(x,0,sigma)+(1-probx)*dnorm(x,mux,sigma),lwd=2.5)
- hist(mu.s,nclass=20,probability=T,xlab='mu')
- hist(prob.s,nclass=20,probability=T,xlab='pi')
- persp(hist2d(mu.s,prob.s,nxbins=17,nybins=17),ylab='pi',xlab='mu',zlab='')
- par(mfrow=c(1,1),bty='n',oma=c(4,4,4,4),pty='s')
- plot(y,pz,xlim=c(-4,6),ylab="estimated class probs",xlab="ordered data y")
# plot all saved MCMC iterations to check convergence
- par(mfrow=c(2,1),bty='n',pty='r',oma=c(1,1,1,1))
- tsplot(mu.s,xlab="MCMC iteration",ylab="mu")
- tsplot(prob.s,xlab="MCMC iteration",ylab="pi"
-
- dev.off()