# Question 3a

log.dens <- function(x1,x2,nu) {             # log of the joint density of x1,x2
  a <- (x1+x2)/(nu+2)
  log(nu) + log(nu+1) - (nu+2) * log(exp(a)+exp(a-x1)+exp(a-x2))
}

nu <- 1                                      # set parameter nu  

####################################################

t <- seq(-5, 5, length.out=200)              # define one side of the grid
mat <- exp(outer(t,t,"log.dens",nu=nu))      # evaluate the density over a grid defined by t
filled.contour(t,t,mat)                      # create the plot
title("Density for nu=1")                    # give it a title

####################################################

n <- 1e4                                     # set the sample size
x <- matrix(nrow=n,ncol=2)                   # create matrix to hold sample
cur.x <- c(0,0)                              # set the initial value
cur.logf <- log.dens(cur.x[1],cur.x[2],nu=nu)# evaluate density at cur.x
n.accepted <- 0                              # counter of the accepted values
for (i in 1:n) {                             # repeat n  times ...
  new.x <- cur.x + sqrt(9) * rnorm(2)        # create proposed value
  new.logf <- log.dens(new.x[1],new.x[2],nu=nu)
                                             # evaluate density at proposed value
  prob.acceptance <- exp(new.logf-cur.logf)  # probability of acceptance
  if (runif(1)<prob.acceptance) {            # if we accept the proposed value ...
    n.accepted <- n.accepted+1               # ... increment the counter of accepted values}
    cur.x <- new.x                           # ... accept the new value
    cur.logf <- new.logf                     # ... retain the density at the accepted value
  }
  x[i,] <- cur.x                             # store current value
}

####################################################

par(mfrow=c(1,3))                            # split the screen into three
plot(x)                                      # plot the sample
title("Joint distribution")                  # add a title
plot(x[,1],type="l")                         # plot the sample path of x1
title("Sample path of X1")                   # add a title
plot(x[,2],type="l")                         # plot the sample path of x2
title("Sample path of X2")                   # add a title

####################################################

mean(x[,1])                                  # Monte Carlo estimate of E(X1)
mean(x[,2])                                  # Monte Carlo estimate of E(X2)

####################################################

n.accepted/n                                 # proportion of accepted values

####################################################

cor.x1 <- cor(x[-1,1],x[-nrow(x),1])         # autocorrelation of x1
ess.x1 <- n * (1-cor.x1) / (1+cor.x1)        # corresponding effective sample size
cor.x2 <- cor(x[-1,2],x[-nrow(x),2])         # autocorrelation of x2
ess.x2 <- n * (1-cor.x2) / (1+cor.x2)        # corresponding effective sample size

####################################################

# Question 3b

sample.cond.dist <- function(x2,nu) {        # function that samples from x1 given x2
  u <- runif(1)                              # generate uniform u
  v <- u^(1/(nu+1))                          # compute intermediate result
  -log((1+exp(-x2))*(1-v)/v)                 # return inverse CDF at u
}

####################################################

x <- matrix(nrow=n,ncol=2)                   # create matrix to hold sample
cur.x <- c(0, 0)                             # set the initial value
for (i in 1:n) {                             # repeat n  times ...
  cur.x[1] <- sample.cond.dist(cur.x[2],nu=nu)
                                             # update x1 by sampling from x1 given x2         
  cur.x[2] <- sample.cond.dist(cur.x[1],nu=nu)
                                             # update x2 by sampling from x2 given x1 
  x[i,] <- cur.x                             # store current value
}

####################################################

# Question 3c

nu <- 0.05                                   # set parameter nu  
t <- seq(-50, 10, length.out=200)            # define one side of the grid
mat <- exp(outer(t,t,"log.dens",nu=nu))      # evaluate the density over a grid defined by t
filled.contour(t,t,mat)                      # create the plot
title("Density for nu=0.05")                 # give it a title

####################################################

library(MASS)                                # load the library MASS (required for mvrnorm)
sigma.proposal <- matrix(c(1000,998,998,1000),nrow=2)
                                             # set covariance of the proposal
x <- matrix(nrow=n,ncol=2)                   # create matrix to hold sample
cur.x <- c(0,0)                              # set the initial value
cur.logf <- log.dens(cur.x[1],cur.x[2],nu=nu)# evaluate density at cur.x
n.accepted <- 0                              # counter of the accepted values
for (i in 1:n) {                             # repeat n  times ...
  new.x <- cur.x + mvrnorm(1,mu=c(0,0),Sigma=sigma.proposal)
                                             # create proposed value with the above covariance
  new.logf <- log.dens(new.x[1],new.x[2],nu=nu)     
                                             # evaluate density at proposed value
  prob.acceptance <- exp(new.logf-cur.logf)  # probability of acceptance
  if (runif(1)<prob.acceptance) {            # if we accept the proposed value ...
    n.accepted <- n.accepted+1               # ... increment the counter of accepted values
    cur.x <- new.x                           # ... accept the new value
    cur.logf <- new.logf                           # ... retain the density at the accepted value
  }
  x[i,] <- cur.x                             # store current value
}

####################################################

# Question 4

load("/home/ioana/public_html/mixgamma.RData")                 

####################################################

log.posterior <- function(lambda, x) {           # evaluates the log-posterior
  sum(log(0.5*dgamma(x,4,lambda[1])+0.5*dgamma(x,5,lambda[2])))+sum(dgamma(lambda,4,0.5,log=TRUE))
}

####################################################

sd.proposal <- 1                             # use large enough a standard deviation of the proposal
n <- 1e4                                     # set the sample size
lambda <- matrix(ncol=2,nrow=n)              # create matrix to hold sample
cur.lambda <- c(2,8)                         # set the initial value
cur.logf <- log.posterior(cur.lambda, x) + sum(log(cur.lambda))
                                             # evaluate log-posterior at cur.lambda
n.accepted <- 0                              # counter of the accepted values
for (i in 1:n) {                             # repeat n  times ...
  new.lambda <- cur.lambda * exp(sd.proposal*rnorm(2))
                                             # create proposed value using a multiplicative perturbation}
  new.logf <- log.posterior(new.lambda, x) + sum(log(new.lambda))
                                             # evaluate log-posterior at proposed value
  prob.acceptance <- exp(new.logf-cur.logf)  # probability of acceptance
  if (runif(1)<prob.acceptance) {            # if we accept the proposed value ...  
    n.accepted <- n.accepted+1               # ... increment the counter of accepted values
    cur.lambda <- new.lambda                 # ... accept the new value
    cur.logf <- new.logf                     # ... retain the density at the accepted value
  }
  lambda[i,] <- cur.lambda                   # store current value
}

####################################################

par(mfrow=c(1,3))                            # split the screen into three
hist(lambda[,1])                             # create a histogram of lambda1
hist(lambda[,2])                             # create a histogram of lambda2
plot(lambda)                                 # create a plot of the joint distribution
title("Joint distribution")                  # add a title

####################################################

# Question 5

log.dens <- function(x1, x2, lambda=5) {     # define the density as R function
  -lambda*(abs(x1^2+x2^2-1))
}

####################################################

lambda <- 5                                  # define lambda
t <- seq(-1.5, 1.5, length.out=200)          # define one side of the grid
mat <- exp(outer(t,t,"log.dens",lambda=lambda))
                                             # evaluate the density over a grid defined by t
filled.contour(t,t,mat)                      # create the plot

####################################################

sd.proposal <- 0.2                           # set standard deviation of the proposal
n <- 1e4                                     # set the sample size
x <- matrix(nrow=n,ncol=2)                   # create matrix to hold sample
cur.x <- c(0,0)                              # set the initial value
cur.logf <- log.dens(cur.x[1],cur.x[2])      # evaluate log-density at cur.x
n.accepted <- 0                              # counter of the accepted values
for (i in 1:n) {                             # repeat n  times ...             
  new.x <- cur.x + sd.proposal * rnorm(2)    # create proposed value
  new.logf <- log.dens(new.x[1],new.x[2])    # evaluate log-density at proposed value}
  prob.acceptance <- exp(new.logf-cur.logf)  # probability of acceptance
  if (runif(1)<prob.acceptance) {            # if we accept the proposed value ...  
    n.accepted <- n.accepted+1               # ... increment the counter of accepted values
    cur.x <- new.x                           # ... accept the new value
    cur.logf <- new.logf                     # ... retain the density at the accepted value
  } 
  x[i,] <- cur.x                             # store current value
}
n.accepted/n                                 # proportion of accepted values
cor(x[-1,1],x[-nrow(x),1])                   # autocorrelation of X1
plot(x,cex=0.2)                              # plot the sample

####################################################

sd.proposal <- 0.1                           # set standard deviation of the proposal
n <- 1e4                                     # set the sample size
x <- matrix(nrow=n,ncol=2)                   # create matrix to hold sample
cur.x <- c(0,0)                              # set the initial value
cur.logf <- log.dens(cur.x[1],cur.x[2])      # evaluate density at cur.x
n.accepted <- 0                              # counter of the accepted values
for (i in 1:n) {                             # repeat n  times ...              
  new.x <- cur.x + sd.proposal * rnorm(2)    # MH MOVE 1: create proposed value
  new.logf <- log.dens(new.x[1],new.x[2])    # MH MOVE 1: evaluate density at proposed value
  prob.acceptance <- exp(new.logf-cur.logf)  # MH MOVE 1: probability of acceptance
  if (runif(1)<prob.acceptance) {            # MH MOVE 1: if we accept the proposed value ...
    n.accepted <- n.accepted+1               # ... increment the counter of accepted values
    cur.x <- new.x                           # ... accept the new value
    cur.logf <- new.logf                     # ... retain the density at the accepted value
  } 
  angle <- 2*pi*runif(1)                     # MH MOVE 2: draw an angle at random
  rotation.matrix <- matrix(c(cos(angle),-sin(angle),sin(angle),cos(angle)),nrow=2)
                                             # MH MOVE 2: Rotation matrix
  cur.x <- rotation.matrix%*%cur.x           # MH MOVE 2: perform the rotation
  x[i,] <- cur.x                             # store current value
}
n.accepted/n                                 # proportion of accepted values
cor(x[-1,1],x[-nrow(x),1])                   # autocorrelation of X1
plot(x,cex=0.2)                              # plot the sample

####################################################

n <- 1e4                                     # set the sample size
u <- runif(n)                                # generate n uniforms
right.tail.prob <- 1/(2-exp(-5))             # compute probability of r greater than 1  
r <- ifelse(runif(n)<right.tail.prob,1-log(u)/5,1+log(exp(-5)+(1-exp(-5))*u)/5)
                                             # sample r (n times)
theta <- 2*pi*runif(n)                       # sample angle theta (n times)
x <- cbind(sqrt(r)*cos(theta),sqrt(r)*sin(theta))
                                             # transform to x1 and x2
plot(x,cex=0.2)                              # plot the results to ensure everything went OK

####################################################

# Question 6c

leaf <- rbind(c(3.28, 3.09, 3.03, 3.03), 
              c(3.52, 3.48, 3.38, 3.38),
              c(2.88, 2.80, 2.81, 2.76),
              c(3.34, 3.38, 3.24, 3.26))
y.total <- mean(leaf)			# compute mean of all observations
y.bar <- apply(leaf, 1, mean)		# compute mean of rows in observation matrix
n <- 1e5				# set sample size
alphas <- matrix(nrow=n, ncol=4)	# create matrix and vectors to hold samples
mu <- numeric(n)
sigma2.alpha <- numeric(n)
sigma2 <- numeric(n)
cur.alpha <- rep(0.01, times=4)		# set starting values
cur.sigma2 <- 0.005
for(i in 1:n) {                         # Gibbs sampler algorithm
  cur.mu <- y.total - mean(cur.alpha) + sqrt(cur.sigma2/16) * rnorm(1)
                                        # update mu
  sigma2.alpha.rate <- 0.5 * sum(cur.alpha^2)
  cur.sigma2.alpha <- 1/rgamma(1, 1, sigma2.alpha.rate)
                                        # update sigma2.alpha
  sigma2.rate <- 0.5 * sum( (leaf - cur.alpha - cur.mu)^2 )
  cur.sigma2 <- 1/rgamma(1, 7, sigma2.rate) 
                                        # update sigma2
  alpha.ratio <- cur.sigma2.alpha / (4 * cur.sigma2.alpha + cur.sigma2)
  cur.alpha <- 4*alpha.ratio*(y.bar-cur.mu) + sqrt(cur.sigma2 * alpha.ratio) * rnorm(4)
                                        # update alphas
  alphas[i,] <- cur.alpha               # Store current values
  mu[i] <- cur.mu
  sigma2[i] <- cur.sigma2
  sigma2.alpha[i] <- cur.sigma2.alpha
}

####################################################

s <- seq(60001,100000, by=2)
par(mfrow=c(1,2))
plot(sigma2[s], type="l", xlab="t", ylab=expression((sigma^2)[t]))
plot(sigma2.alpha[s], type="l", xlab="t", ylab=expression((sigma[alpha]^2)[t]))

par(mfrow=c(1,2))
hist(sigma2[s], xlab=expression(sigma^2), main="",freq=FALSE)
hist(sigma2.alpha[s], xlab=expression(sigma[alpha]^2), main="",freq=FALSE)

####################################################

s1 <- seq(60003, 100000, by=2)
s2 <- seq(60001, 99997, by=2)
cor(sigma2[s1], sigma2[s2])
cor(sigma2.alpha[s1], sigma2.alpha[s2])
s.all <- 60000:100000
median(sigma2[s.all])
median(sigma2.alpha[s.all])