#beta posterior
N = 1000
theta = numeric(N)
theta[1] = runif(1)
for(i in 1:(N-1)) {
  theta.star = runif(1, theta[i]-.2, theta[i]+.2)
  if(theta.star >= 0 & theta.star <= 1) {
    alpha = theta.star^6 * (1 - theta.star)^4 / (theta[i]^6 * (1 - theta[i])^4)
    if(runif(1) < alpha) theta[i+1] = theta.star
    else theta[i+1] = theta[i]
  }
  else theta[i+1] = theta[i]
}

hist(theta)
x=seq(0,1,by=.01)
lines(x,dbeta(x,7,5)*N*.05)

#trace plot
plot(theta,pch=".")
mean(theta)
var(theta)

N = 10000

#track acceptance probability
N = 1000
theta = numeric(N)
theta[1] = runif(1)
accept = 0
for(i in 1:(N-1)) {
  theta.star = runif(1, theta[i]-.2, theta[i]+.2)
  if(theta.star >= 0 & theta.star <= 1) {
    alpha = theta.star^6 * (1 - theta.star)^4 / (theta[i]^6 * (1 - theta[i])^4)
    if(runif(1) < alpha) {
      accept = accept + 1
      theta[i+1] = theta.star
    }
    else theta[i+1] = theta[i]
  }
  else theta[i+1] = theta[i]
}
accept/N

#change proposal function
theta = numeric(N)
theta[1] = runif(1)
accept = 0
for(i in 1:(N-1)) {
  theta.star = runif(1)
  alpha = theta.star^6 * (1 - theta.star)^4 / (theta[i]^6 * (1 - theta[i])^4)
  if(runif(1) < alpha) {
    accept = accept + 1
    theta[i+1] = theta.star
  }
  else theta[i+1] = theta[i]
}
accept/N
hist(theta)
plot(theta,pch=".")

#acceptance rate high, mixing slower
theta = numeric(N)
theta[1] = runif(1)
accept = 0
for(i in 1:(N-1)) {
  theta.star = runif(1, theta[i]-.1, theta[i]+.1)
  if(theta.star >= 0 & theta.star <= 1) {
    alpha = theta.star^6 * (1 - theta.star)^4 / (theta[i]^6 * (1 - theta[i])^4)
    if(runif(1) < alpha) {
      accept = accept + 1
      theta[i+1] = theta.star
    }
    else theta[i+1] = theta[i]
  }
  else theta[i+1] = theta[i]
}
accept/N
plot(theta,pch=".")


#normal likelihood with Cauchy prior
N = 1000
y.bar = 5.4
n = 20
mu = numeric(N)
mu[1] = rnorm(1)
for(i in 1:(N-1)) {
  mu.star = rnorm(1,mu[i])
  alpha = exp( (n/2) * ((mu[i] - y.bar)^2 - (mu.star - y.bar)^2) ) *
     (1 + mu[i]^2) /(1 + mu.star^2)
  if(runif(1) < alpha) mu[i+1] = mu.star
  else mu[i+1] = mu[i]
}
hist(mu)
plot(mu,pch=".")

#remove burn-in and track acceptance probability
N.iter = 1000
N.burnin = 50
y.bar = 5.4
n = 20
mu = numeric(N.iter+N.burnin)
mu[1] = rnorm(1)
accept=0
for(i in 1:(N.iter+N.burnin-1)) {
  mu.star = rnorm(1,mu[i])
  alpha = exp( (n/2) * ((mu[i] - y.bar)^2 - (mu.star - y.bar)^2) ) *
     (1 + mu[i]^2) /(1 + mu.star^2)
  if(runif(1) < alpha) {
    mu[i+1] = mu.star
    accept = accept + 1
  }
  else mu[i+1] = mu[i]
}
mu = mu[(N.burnin+1):(N.iter+N.burnin)]
hist(mu)
plot(mu,pch=".")
accept/(N.iter+N.burnin)
mean(mu)
var(mu)


#normal likelihood with mean and variance unknown
N.iter = 1000
N.burnin = 100
y = rnorm(20, mean=5.4, sd=2.1)
y.bar = mean(y)
n = 20
mu = numeric(N.iter+N.burnin)
sigma2 = numeric(N.iter+N.burnin)
mu[1] = rnorm(1)
sigma2[1] = rexp(1)
accept.mu = 0
accept.sigma2 = 0
for(i in 1:(N.iter+N.burnin-1)) {
  #sample mu using a normal proposal
  mu.star = rnorm(1,mu[i],sd=sqrt(sigma2[i]/n))
  alpha = exp((n/(2*sigma2[i]^2)) * ((mu[i] - y.bar)^2 - (mu.star - y.bar)^2))*
     (1 + mu[i]^2) /(1 + mu.star^2)
  if(runif(1) < alpha) {
    mu[i+1] = mu.star
    accept.mu = accept.mu + 1
  }
  else mu[i+1] = mu[i]
  #sample sigma^2 using an asymmetric uniform proposal
  sigma2.star = runif(1,sigma2[i]/2,sigma2[i]*2)
  ysum = sum((y-mu[i+1])^2)  #could also compute this as t(y-mu[i])%*%(y-mu[i])
  alpha = exp(ysum/(2*sigma2[i])+sigma2[i]-ysum/(2*sigma2.star)-sigma2.star)*
    sigma2[i]^(n/2-1)/sigma2.star^(n/2-1)
  if(runif(1) < alpha) {
    sigma2[i+1] = sigma2.star
    accept.sigma2 = accept.sigma2 + 1
  }
  else sigma2[i+1] = sigma2[i]
}
mu = mu[(N.burnin+1):(N.iter+N.burnin)]
sigma2 = sigma2[(N.burnin+1):(N.iter+N.burnin)]
hist(mu)
plot(mu,pch=".")
accept.mu/(N.iter+N.burnin)
mean(mu)
var(mu)

hist(sigma2)
plot(sigma2,pch=".")
accept.sigma2/(N.iter+N.burnin)
mean(sigma2)
var(sigma2)