source("lab6.R")
##########################################################
##Examples
##########################################################
## Independence sampler - poor proposal distribution
##########################################################
f=independence.sampler.markov.chain(
  function(x) -100*(x< -1 | x>1), #log.pi - this is basically the uniform[-1,1] distribution
  function() rnorm(1,sd=0.1), #Normal(0,0.1^2) proposal - poor choice!
  function(x) log(dnorm(x,sd=0.1)),
  iterations=10000,
  burn.in=1000,
  thin=1)
chain=f(0)
raftery.diag(chain)
summary(chain)
plot(chain)
geweke.diag(chain)
geweke.plot(chain)
run.gelman.diagnostic(runif(10,-1,1),f)


##########################################################
## Independence sampler - better proposal
##########################################################
f=independence.sampler.markov.chain(
  function(x) -100*(x< -1 | x>1), #log.pi - this is basically the uniform[-1,1] distribution
  function() rnorm(1,sd=1), 
  function(x) log(dnorm(x,sd=1)),
  iterations=10000,
  burn.in=1000,
  thin=1)
chain=f(0)
raftery.diag(chain)
summary(chain)
plot(chain)
geweke.diag(chain)
geweke.plot(chain)
run.gelman.diagnostic(runif(10,-1,1),f)

##########################################################
## Independence samples - better proposal
##########################################################
f=independence.sampler.markov.chain( ##2d multivariate normal example
  function(x) -x[1]^2-2*x[1]*x[2]-x[2]^2,
  function() rnorm(2,sd=3),
  function(x) sum(log(dnorm(x,sd=3))))
chain=f(c(0,0))
plot(as.matrix(chain),asp=1)
raftery.diag(chain)
summary(chain)
plot(chain)
plot(as.matrix(chain),asp=1)
geweke.diag(chain)
geweke.plot(chain)
run.gelman.diagnostic(lapply(1:10,function(x) rnorm(8,sd=3)),f)


##########################################################
## Random walk sampler
##########################################################
f=random.walk.sampler.markov.chain(
  function(x) -x^2/(2*10^2), #Normal(0,10^2)
  1,
  iterations=100000)
chain=f(0)
plot(as.matrix(chain),asp=1)
raftery.diag(chain)
summary(chain)
plot(chain)
plot(as.matrix(chain),asp=1)
geweke.diag(chain)
geweke.plot(chain)
run.gelman.diagnostic(lapply(1:10,function(x) rnorm(1,sd=20)),f)

##########################################################
## Random walk sampler
##########################################################
f=random.walk.sampler.markov.chain( #2d multivariate normal example
  function(x) exp(-x[1]^2-2*x[1]*x[2]-x[2]^2),
  3)
chain=f(c(0,0))
plot(as.matrix(chain),asp=1)
raftery.diag(chain)
summary(chain)
plot(chain)
plot(as.matrix(chain),asp=1)
geweke.diag(chain)
geweke.plot(chain)
run.gelman.diagnostic(lapply(1:10,function(x) rnorm(2,sd=3)),f)

##########################################################
## Random walk sampler
##########################################################
f=random.walk.sampler.markov.chain(
  function(x) 100*(abs(sum(x^2)-9)<1), #An annulus, centered on the origin, radius 3 approximately.
  0.3,  
  iterations=10000,
  burn.in=1000,
  thin=1)
chain=f(c(3,0))
plot(as.matrix(chain),asp=1)
raftery.diag(chain)
summary(chain)
plot(chain)
geweke.diag(chain)
geweke.plot(chain)
run.gelman.diagnostic(lapply(1:10,function(x) rnorm(2,sd=3)),f)

##########################################################
## Gibbs sampler
##########################################################
## Pick 4 points in [0,1] uniformly at random, but 
## conditioned to be in increasing order and at distance at 
## least 0.1 from each other.

resample.coordinate <- function(i) {
  function(x) { #Sample x[i] conditional on the rest
    if (i==1)
      min=0
    else
      min=x[i-1]+0.1
    if (i==4)
      max=1
    else
      max=x[i+1]-0.1
    runif(1,min,max)
  }
}

f=gibbs.sampler.markov.chain(
  list(resample.coordinate(1),
       resample.coordinate(2),
       resample.coordinate(3),
       resample.coordinate(4)))
chain=f(c(0,0.3,0.6,1))
summary(chain)
plot(chain)

