###################################################
### chunk number 1: 
###################################################
options(keep.source=TRUE,continue=" ",prompt=" ")


###################################################
### chunk number 2: 
###################################################
  mystery.fun <- function (theta) {
    20*dnorm(x=theta,mean=-1,sd=5)+ifelse(theta>-1,6*dgamma(x=theta+1,shape=2,scale=0.5),dgamma(x=-theta,shape=5,scale=0.2))+2*dgamma(x=3-theta,shape=4,scale=0.25)
  }


###################################################
### chunk number 3: 
###################################################
require(coda)

mcmc.mh <- function (theta, fun, niter, thin = 1) {
  chain <- matrix(
                  nrow=as.integer(niter),
                  ncol=length(theta),
                  dimnames=list(NULL,names(theta))
                  )
  f <- fun(theta)
  chain[1,] <- theta
  k <- 1                                # step in the chain
  nsave <- 1                            # number of values saved so far
  nprop <- 0                            # number of proposals made
  nacpt <- 0                            # number of proposals accepted
  while (nsave < niter) {
    theta.star <- runif(n=1,min=theta-10,max=theta+10) # propose a new parameter
    f.star <- fun(theta.star)
    nprop <- nprop+1
    alpha <- f.star/f  # Metropolis' rule (since proposal is symmetric)
    if ((alpha>=1)||(runif(1)<alpha)) {
      theta <- theta.star
      f <- f.star
      nacpt <- nacpt+1
    }
    k <- k+1
    if (k%%thin==0) {
      nsave <- nsave+1
      chain[nsave,] <- theta
    }
  }
  cat("acceptance fraction = ",nacpt/nprop,"\n")
  mcmc(data=chain,thin=thin)            # store the result in a CODA object
}


###################################################
### chunk number 4: 
###################################################
theta <- mcmc.mh(theta=10,fun=mystery.fun,niter=100)
plot(theta)


###################################################
### chunk number 5: 
###################################################
theta <- mcmc.mh(theta=10,fun=mystery.fun,niter=10000)


###################################################
### chunk number 6: 
###################################################
effectiveSize(theta)
raftery.diag(theta)


###################################################
### chunk number 7: 
###################################################
theta <- mcmc.mh(theta=10,fun=mystery.fun,niter=1000)  ## burn-in: discard these
theta <- mcmc.mh(theta=tail(theta,1),fun=mystery.fun,niter=16000,thin=5)


###################################################
### chunk number 8: 
###################################################
effectiveSize(theta)


###################################################
### chunk number 9: 
###################################################
pacf(theta)


###################################################
### chunk number 10: 
###################################################
theta <- mcmc.mh(theta=tail(theta,1),fun=mystery.fun,niter=10000,thin=10)


###################################################
### chunk number 11: 
###################################################
geweke.diag(theta)
pnorm(q=-abs(geweke.diag(theta)$z))
effectiveSize(theta)
heidel.diag(theta)


###################################################
### chunk number 12: 
###################################################
theta <- mcmc.mh(theta=tail(theta,1),fun=mystery.fun,niter=25000,thin=20)


###################################################
### chunk number 13: 
###################################################
heidel.diag(theta)
geweke.diag(theta)
pnorm(q=-abs(geweke.diag(theta)$z))
effectiveSize(theta)


###################################################
### chunk number 14: 
###################################################
pacf(theta)


###################################################
### chunk number 15: 
###################################################
hist(theta,breaks=250,freq=F)
z <- integrate(mystery.fun,lower=-20,upper=20)$value
curve(1/z*mystery.fun(theta=x),from=-20,to=20,col='red',add=T,n=1001)


###################################################
### chunk number 16: 
###################################################
worm.data <- log(1+round(ifelse(runif(1000)<0.2,rlnorm(n=1000,meanlog=log(1000),sdlog=0.5),rlnorm(n=1000,meanlog=log(50),sdlog=1))))


###################################################
### chunk number 17: 
###################################################
hist(worm.data,breaks=50,freq=T,xlab="log worm load")


###################################################
### chunk number 18: 
###################################################
gibbs.mixture <- function (theta, X, niter, thin = 1) {
  ## use diffuse priors:
  ## mu ~ normal(m,1/sqrt(q)), m = 0, q = 0.0001
  ## p ~ uniform(0,1) = beta(a,b), a = b = 1
  prior.m <- c(mu.1=0,mu.2=0)           # mean of normal priors
  prior.q <- c(q1=0.0001,q2=0.0001)     # precision of normal priors
  prior.ab <- c(a=1,b=1)                # beta dist. shape parameters
  true.sigma <- c(sigma1=0.5,sigma2=1)  # normal SDs assumed known
  true.q <- 1/true.sigma^2              # precisions
  chain <- matrix(
                  nrow=as.integer(niter),
                  ncol=length(theta),
                  dimnames=list(NULL,c("p","mu.1","mu.2"))
                  )
  p <- theta[1]
  mu <- theta[-1]
  chain[1,] <- c(p,mu)
  k <- 1
  nsave <- 1
  while (nsave < niter) {
    ## first we choose Z | p,mu
    f1 <- dnorm(x=X,mean=mu[1],sd=true.sigma[1]) # Prob[X | Z==1,mu]
    f2 <- dnorm(x=X,mean=mu[2],sd=true.sigma[2]) # Prob[X | Z==2,mu]
    prob1 <- p*f1/(p*f1+(1-p)*f2)                # Prob[Z=1 | mu,p,X]
    Z <- ifelse(runif(n=length(X))<prob1,1,2)
    ## now choose mu | p,Z
    n <- c(sum(Z==1),sum(Z==2))
    s <- c(sum(X[Z==1]),sum(X[Z==2]))
    post.q <- prior.q+n*true.q
    post.m <- (prior.q*prior.m+true.q*s)/post.q
    mu <- rnorm(n=2,mean=post.m,sd=1/sqrt(post.q))
    ## now choose p | Z
    post.ab <- prior.ab+n
    p <- rbeta(n=1,shape1=post.ab[1],shape2=post.ab[2])
    k <- k+1
    if (k%%thin==0) {
      nsave <- nsave+1
      chain[nsave,] <- c(p,mu)
    }
  }
  mcmc(data=chain,thin=thin)            # store the result in a CODA object
}


###################################################
### chunk number 19: 
###################################################
theta <- gibbs.mixture(theta=c(p=0.5,mu=c(0,0)),X=worm.data,niter=20)


###################################################
### chunk number 20: 
###################################################
plot(theta)


###################################################
### chunk number 21: 
###################################################
theta <- gibbs.mixture(theta=c(p=0.5,mu=c(7,4)),X=worm.data,niter=4000)


###################################################
### chunk number 22: 
###################################################
effectiveSize(theta)
raftery.diag(theta)


###################################################
### chunk number 23: 
###################################################
theta <- gibbs.mixture(theta=theta[nrow(theta),],X=worm.data,niter=5000,thin=4)


###################################################
### chunk number 24: 
###################################################
heidel.diag(theta)
geweke.diag(theta)
pnorm(q=-abs(geweke.diag(theta)$z))


###################################################
### chunk number 25: 
###################################################
plot(theta)


###################################################
### chunk number 26: 
###################################################
summary(theta)


###################################################
### chunk number 27: 
###################################################
theta <- gibbs.mixture(theta=c(p=0.5,mu=c(4,7)),X=worm.data,niter=4000)
theta <- gibbs.mixture(theta=theta[nrow(theta),],X=worm.data,niter=10000,thin=4)


###################################################
### chunk number 28: 
###################################################
plot(theta)


###################################################
### chunk number 29: 
###################################################
summary(theta)