# 1. Accept/reject for target=Normal and proposal: G=2-sided exponential.

# Explaining the underlying principle
par(mfrow=c(2,1))
z <- seq(-5,5,0.01)
f <- dnorm(z) # target density
g <- 0.5*exp(-abs(z)) # proposal density
acc <- sqrt(pi*0.5)*exp(-0.5)  # acceptance probability = 1/M
plot(z,g,type="l",ylab="density",main=
     "normal and 2-sided exponential density") 
abline(h=0)
lines(z,f,col=3) 
lines(z,acc*f,col=2) # f/M
plot(z,acc*f/g,type="l",ylab="f(x)/(g(x)M)",main=
     "acceptance probabilities")

# Implementing 
par(mfrow=c(2,1))
x <- (1-2*(runif(10000) < 0.5))*log(runif(10000)) # Sample from G
qqnorm(x) # G has heavier tails than F
a <- exp(-0.5*(abs(x)-1)^2) # vector of acceptance probabilities
y <- x[runif(10000)<a] # Accept/reject step
qqnorm(y) # Check for normality 
length(y) # Number of accepted values

br <- c(-12,seq(-4,4,0.2),12) # histogram boundaries 
hg <- hist(x,breaks=br,plot=F)$counts # Frequencies of proposed values
hf <- hist(y,breaks=br,plot=F)$counts # Frequencies of accepted values
par(mfrow=c(2,1))
barplot(height=rbind(hf,hg-hf),main=
        "Histogram of proposed and accepted values")
plot(z,g,type="l",ylab="Dichte",main=
     "normal and 2-sided exponential density") 
abline(h=0)
lines(z,f,col=3) 
lines(z,acc*f,col=2)

# 2 Importance sampling for the same pair target/proposal
x <- (1-2*(runif(10000) < 0.5))*log(runif(10000)) #Sample from G
gew <- exp(-0.5*(abs(x)-1)^2)
gew <- gew/sum(gew) #Weights for importance sampling
ESS <- 1/sum(gew^2) # Equivalent sample size
plot(x,gew,type="h",ylab="weights")
plot(sort(gew),ylab="sorted weights",type="l")
abline(h=0.0001,col=2) #shows non-uniformity of weights.

par(mfrow=c(1,1)) #Histogram of the weighted sample
kmax <- ceiling(10*(max(x)-min(x))) #number of classes with length 0.1
c0 <- min(x) - 0.05 
c1 <- c0 + kmax*0.1
z <- ceiling((x-c0)*10) # assigning observations to intervals
h <- rep(0,kmax)
for (k in (1:kmax)) h[k] <- sum(gew[z==k]) # counting frequencies
plot(c0+((1:kmax)-0.5)*0.1,h*10,type="h",xlab="x",ylab="density" )
xx <- seq(-3,3,0.01)
lines(xx,dnorm(xx),col=2)


# 3. i.i.d normals with non-conjugate prior

#generating data
x <- rnorm(4) 
xb <- mean(x)
n <- length(x)
xs <- var(x)*(n-1)/n

#Choices  of hyperparameter
#a). not informative
xi <- 0
ka2 <- 100*xs
gam <- 1
lam <- 1
#b) informative for theta, in agreement with data
xi <- xb
ka2 <- xs/n
gam <- 1
lam <- 1
#c) informative for theta, conflicting the data
xi <- xb+3*sqrt(xs/n)
ka2 <- xs/n
gam <- 1
lam <- 1

f.post.contour <- function(xi,ka2,gam,lam,n,xb,xs)
{
  ## Purpose: Plotting contours of the posterior. The main task is to define
  ##          a suitable domain where the posterior is concentrated 
  ## -------------------------------------------------------------------------
  ## Arguments: xi,ka2,al,lam=hyperparameters, xb,xs=mle
  ##            n=sample size (n=0 gives the prior !)
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date:  1 Nov 1999, 18:26
  s2 <- lam/gam #gam/lam is the prior mean of 1/sigma^2
  a <- s2/(s2+n*ka2)
  mu <- a*xi+(1-a)*xb # posterior mean of theta if sigma^2=s2
  th <- seq(mu-3*sqrt(a*ka2),mu+3*sqrt(a*ka2),length=32)
                                        # defines the range for theta
  lam <- lam+0.5*n*xs 
  gam <- gam+0.5*n
     #if theta=xb, posterior for 1/sigma^2 is Gamma(lam,gam)
  si <- sqrt(lam)*seq(1/sqrt(qgamma(0.99,gam)),1/sqrt(qgamma(0.01,gam)),
                      length=32) # defines the range for sigma 
  z <- 0.5*outer((th-xi)^2,rep(1,32))/ka2+(2*gam+1)*outer(rep(1,32),log(si))+
    outer(lam+0.5*n*(th-xb)^2,si^2,FUN="/")
           #log posterior on the lattice of values for theta and sigma 
  contour(th,si,-z,nlevels=20,xlab="mu",ylab="sigma")
  }

f.lik.contour <- function(n,xb,xs)
{
  ## Purpose: Plotting the log likelihood of normal observations
  ## ----------------------------------------------------------------------
  ## Arguments: xb,xs=mle, n=sample size
  ## ----------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 16 Nov 2003, 23:05
  th <- seq(xb-3*sqrt(xs/n),xb+3*sqrt(xs/n),length=32)
                                        # range for theta
  si <- sqrt(n*xs)*seq(1/sqrt(qchisq(0.99,n-1)),1/sqrt(qchisq(0.01,n-1)),
                      length=32)# range for sigma
  z <- n*outer(rep(1,32),log(si))+ 0.5*n*outer(xs+(th-xb)^2,si^2,FUN="/")
                                 #computing log likelihood
  contour(th,si,-z,nlevels=20,xlab="mu",ylab="sigma")
}

#prior density, likelihood und posterior density 
par (mfrow=c(2,2))
f.post.contour(xi,ka2,gam,lam,0,xb,xs)
title(main="log prior density")
f.lik.contour(n,xb,xs)
title(main="log likelihood")
points(xb,sqrt(xs),pch=15)
f.post.contour(xi,ka2,gam,lam,n,xb,xs)
points(xb,sqrt(xs),pch=15)
title(main="log posterior density")

# 4. Gibbs sampler for i.i.d. normals with non-conjugate prior

f.gibbs <- function(xi,ka2,gam,lam,n,xb,xs,nrep=1000)
{
  ## Purpose: Implementing the Gibbs-sampler for the posterior in a normal
  ## model with independent normal-inverse Gamma prior
  ## -------------------------------------------------------------------------
  ## Arguments: xi,ka2,gam,lam=hyperparameters, xb,xs=sufficient statistic
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date:  1 Nov 1999, 15:50
  th <- rep(xb,nrep)
  si <- rep(xs,nrep)
  for (k in (2:nrep)){
    a <- si[k-1]/(si[k-1]+n*ka2)
    th[k] <- rnorm(1,a*xi+(1-a)*xb,sqrt(a*ka2))
    si[k] <- (lam+0.5*n*(xs+(xb-th[k])^2))/rgamma(1,0.5*n+gam)
  }
  list(th=th,si=sqrt(si))
}

# First few iterations of the Gibbs sampler
par(mfrow=c(1,1))
gibbs <- f.gibbs(xi,ka2,gam,lam,n,xb,xs,nrep=50)
f.post.contour(xi,ka2,gam,lam,n,xb,xs)
points(gibbs$th[1],gibbs$si[1],col=4,pch=16)
i <- 0
i <- i+1 # iterate from here to see more than one step
points(gibbs$th[i+1],gibbs$si[i],col=4,pch=16)
arrows(gibbs$th[i],gibbs$si[i],gibbs$th[i+1],gibbs$si[i],length=0.01,col=4)
points(gibbs$th[i+1],gibbs$si[i+1],col=4,pch=16)
arrows(gibbs$th[i+1],gibbs$si[i],gibbs$th[i+1],gibbs$si[i+1],length=0.1,col=4)

# Large sample
gibbs <- f.gibbs(xi,ka2,gam,lam,n,xb,xs,nrep=1000)
f.post.contour(xi,ka2,gam,lam,n,xb,xs)
title(main="Gibbs-Sampler and log posterior density")
points(gibbs$th,gibbs$si)
par(mfrow=c(2,1))
plot(gibbs$th[1:500],type="l",ylab="mu",xlab="iteration")
             #trace plots (time series of each component)
plot(gibbs$si[1:500],type="l",ylab="sigma",xlab="iteration")
acf(gibbs$th,lag.max=100)  #Autocorrelations 
acf(gibbs$si,lag.max=100)
par(mfrow=c(1,1))

# 5. Metropolis for i.i.d normals with non-conjugate prior

f.metrop <- function(x0,sigma,nrep,r.target, ...)
{
  ## Purpose: Metropolis algorithm with Gaussian random walk proposals
  ## -------------------------------------------------------------------------
  ## Arguments: x0 = (vector of) starting values, sigma = vector of
  ##            standard deviations for the random walk, nrep=number of
  ##            replicates, r.target = ratio of target densities
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 22 Jan 2002, 19:21
  d <- length(x0)
  x <- matrix(rnorm(d*nrep),nrow=nrep,ncol=d)
  y <- matrix(0,nrow=nrep-1,ncol=d)
  x <- x%*%diag(sigma,nrow=d)
  x[1,] <- x0
  nrej <- 0
  for (i in (2:nrep)) {
    x[i,] <- x[i-1,]+x[i,]
    y[i-1,] <- x[i,]
    if (runif(1) > r.target(x[i,],x[i-1,],...)) {
      x[i,] <- x[i-1,]
      nrej <- nrej+1
    }
  }
  list(sample=x,nrej=nrej,prop=y)
}

r.post <- function(x1,x2,xi,ka2,gam,lam,n,xb,xs)
{
  ## Purpose: Computing ratio of posteriori densities at x1 and x2
  ##          for an i.i.d. normal model
  ## ----------------------------------------------------------------------
  ## Arguments: x=(mu,log(sigma)) = Parameters of the normal distribution
  ##            xi,ka2,gam,lam = Parameters of the prior
  ##            n=number of data, xb=sample mean, xs=sample variance
  ## ----------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 24 Jan 2004, 17:25
  quot <- 0.5*((x2[1]-xi)^2-(x1[1]-xi)^2)/ka2 + (2*gam+n)*(x2[2]-x1[2]) +
    (lam+0.5*n*(xs+(x2[1]-xb)^2))*exp(-2*x2[2]) -
    (lam+0.5*n*(xs+(x1[1]-xb)^2))*exp(-2*x1[2])
  exp(quot)
}

# Few iterations (with log scale for sigma) to explain the idea
f.post.contour(xi,ka2,gam,lam,n,xb,xs)
points(xb,sqrt(xs),pch=15,col=2)
title(main="log posterior density")
metrop <- f.metrop(x0=c(xb,0.5*log(xs)),sigma=c(0.2,0.1),nrep=21,
                   r.target=r.post,xi,ka2,gam,lam,n,xb,xs)
    # vary values of sigma to see how it influences the results
sample <- metrop$sample
sample[,2] <- exp(sample[,2])
prop <- metrop$prop
prop[,2] <- exp(prop[,2])
points(prop[,1],prop[,2],col=3) # proposed points in green
points(sample[,1],sample[,2],col=2) #accepted points in red
arrows(sample[-21,1],sample[-21,2],prop[,1],prop[,2],length=0.1,col=4)
#Result of many iterations
f.post.contour(xi,ka2,gam,lam,n,xb,xs)
points(xb,sqrt(xs),pch=15,col=2)
title(main="log posterior density")
metrop <- f.metrop(x0=c(xb,0.5*log(xs)),sigma=c(0.6,0.3),nrep=5000,
                   r.target=r.post,xi,ka2,gam,lam,n,xb,xs)
metrop$nrej  # number of rejections
param <- metrop$sample
points(param[,1]+0.04*(runif(5000)-0.5),
       exp(param[,2])+0.04*(runif(5000)-0.5),cex=0.5)
                      # overlapping points made visible by jittering
par(mfrow=c(2,1))
plot(param[1:500,1],type="l",xlab="iteration",ylab="mu")
           #trace plots (time series of each component)
plot(param[1:500,2],type="l",xlab="iteration",ylab="sigma")
acf(param[,1],lag.max=100)  #Autocorrelations 
acf(param[,2],lag.max=100)
par(mfrow=c(1,1))


# 5. Metropolis for a Gaussian mixture

f.contour.mixture <- function(p,mu,scale)
{
  ## Purpose: Plotting the density of a mixture of 2 bivariate normals
  ##          The first normal is the standard normal. 
  ## -------------------------------------------------------------------------
  ## Arguments: p = mixture weight, mu=location of second normal,
  ##            scale = scale of second normal (scalar) 
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 22 Jan 2002, 22:07
  x <- seq(-2.5,4,length=32)
  z <- p*exp(-0.5*outer(x^2,x^2,FUN="+"))+
    (1-p)*exp(-0.5*outer((x-mu[1])^2,(x-mu[2])^2,FUN="+")/scale^2)
  contour(x,x,z,nlevels=20)}


r.mixture <- function(x1,x2,p,mu,scale)
{
  ## Purpose: Computing the ratio of a mixture density at two places x1, x2
  ## -------------------------------------------------------------------------
  ## Arguments: p, mu scale as in the previous function
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 22 Jan 2002, 21:25
  r <- p*exp(-0.5*sum(x1^2)) + (1-p)*exp(-0.5*sum((x1-mu)^2)/scale^2)
  r <- r/(p*exp(-0.5*sum(x2^2)) + (1-p)*exp(-0.5*sum((x2-mu)^2)/scale^2))
  r
}



mu <- c(3,3)
scale <- 0.5
#target density and generated sample
par(mfrow=c(1,1))
f.contour.mixture(0.7,mu,scale) 
erg <- f.metrop(x0=c(0,0),sigma=c(1,1),5000,r.mixture,0.7,mu,scale)
                                        # vary argument sigma!
erg$nrej # number of rejected proposals
x <- erg$sample
points(x[,1],x[,2])
# trace plots and autocorrelations 
par(mfrow=c(2,1))
plot(x[1:500,1],type="l")
plot(x[1:500,2],type="l")
acf(x[,1],lag.max=100)
acf(x[,2],lag.max=100)
par(mfrow=c(1,1))