# TEST3.R - Test of short cut Metropolis on the funnel distribution.

source("../short-met.r")


# Description of the distribution.

rdim <- 9		# Number of radial dimensions
sdv <- 3		# Standar deviation of the log variance


# Numbers of iterations & seed.

seed <- 1

cnt <- 5000
cntsc.ra <- 2.1*cnt


# Compute radial distance.

rdist <- function (x)
{ 
  sqrt(sum(x[2:(rdim+1)]^2))
}


# Log probability density function. 

fun.lpr <- function (x)
{
  dnorm(x[1],0,sdv,log=TRUE) + 
    sum(dnorm(x[2:(rdim+1)],0,sqrt(exp(x[1])),log=TRUE))
}


# Proposal generation.

fun.pf <- function ()
{ 
  rnorm(rdim+1,0,1)
}


# Initial state.

initial.x <- c(0,rep(1,rdim))


# Standard Metropolis, w=0.03.

set.seed(seed)
lpr.evals <- 0
system.time (res.stand.0.03 <- 
   standard.final (initial.x, fun.lpr, fun.pf, 0.03, 
                   rep(c(1000,1000,1000,1000),cnt)) )
R.stand.0.03 <- res.stand.0.03$states
rej.stand.0.03 <- res.stand.0.03$rejected
evals.stand.0.03 <- lpr.evals


# Standard Metropolis, w=0.15.

set.seed(seed)
lpr.evals <- 0
system.time (res.stand.0.15 <- 
   standard.final (initial.x, fun.lpr, fun.pf, 0.15, 
                   rep(c(1000,1000,1000,1000),cnt)) )
R.stand.0.15 <- res.stand.0.15$states
rej.stand.0.15 <- res.stand.0.15$rejected
evals.stand.0.15 <- lpr.evals


# Standard Metropolis, w=0.75.

set.seed(seed)
lpr.evals <- 0
system.time (res.stand.0.75 <- 
   standard.final (initial.x, fun.lpr, fun.pf, 0.75,
                   rep(c(1000,1000,1000,1000),cnt)) )
R.stand.0.75 <- res.stand.0.75$states
rej.stand.0.75 <- res.stand.0.75$rejected
evals.stand.0.75 <- lpr.evals


# Standard Metropolis, w=3.75.

set.seed(seed)
lpr.evals <- 0
system.time (res.stand.3.75 <- 
   standard.final (initial.x, fun.lpr, fun.pf, 3.75, 
                   rep(c(1000,1000,1000,1000),cnt)) )
R.stand.3.75 <- res.stand.3.75$states
rej.stand.3.75 <- res.stand.3.75$rejected
evals.stand.3.75 <- lpr.evals


# Standard Metropolis, with w=0.03 for 1000 updates, w=0.15 for 1000 updates,
# w=0.75 for 1000 updates, and w=3.75 for 1000 updates.

set.seed(seed)
lpr.evals <- 0
system.time (res.stand.all <- 
  standard.final (initial.x, fun.lpr, fun.pf, c(0.03,0.15,0.75,3.75), 
                  rep(c(1000,1000,1000,1000),cnt)) )
R.stand.all <- res.stand.all$states
rej.stand.all <- res.stand.all$rejected
evals.stand.all <- lpr.evals


# Short cut Metropolis, with w=0.03 for 1000 updates, w=0.15 for 1000 updates,
# w=0.75 for 1000 updates, and w=3.75 for 1000 updates, reversing on all 
# rejections or on less than 3 rejections in a group of 40 updates.

set.seed(seed)
lpr.evals <- 0
system.time (res.short.ra <- 
  short.cut.final (initial.x, fun.lpr, fun.pf, c(0.03,0.15,0.75,3.75), 
                   c(40,40,40,40), rep(c(25,25,25,25),cntsc.ra), 
                   c(3,3,3,0), c(40,39,39,39)) )
R.short.ra <- res.short.ra$states
rej.short.ra <- res.short.ra$rejected
copied.short.ra <- res.short.ra$copied
index.short.ra <- res.short.ra$index
step.short.ra <- res.short.ra$stepsize
evals.short.ra <- lpr.evals


# Compute autocorrelation times.

ac.time <- function (x, lag.max=50)
{
  2 * sum(acf(as.numeric(x),plot=F,lag.max=lag.max)$acf) - 1
}

ac.times.m <- c(ac.time(R.stand.0.03[,1],lag.max=1000),
                ac.time(R.stand.0.15[,1],lag.max=100),
                ac.time(R.stand.0.75[,1],lag.max=1000),
                ac.time(R.stand.3.75[,1],lag.max=1000),
                ac.time(R.stand.all[,1]),
                ac.time(R.short.ra[,1]))


# Compute and print table of results.

results <- data.frame (list (
  method  = c("stand.0.03", "stand.0.15", "stand.0.75", "stand.3.75", 
              "stand.all", "short.ra"),
  evals   = c(evals.stand.0.03, evals.stand.0.15, 
              evals.stand.0.75, evals.stand.3.75,
              evals.stand.all, evals.short.ra), 
  rej     = round(c(mean(rej.stand.0.03), mean(rej.stand.0.15), 
                    mean(rej.stand.0.75), mean(rej.stand.3.75), 
                    mean(rej.stand.all), mean(rej.short.ra)),3),
  mean    = round(c(mean(R.stand.0.03[,1]), mean(R.stand.0.15[,1]), 
                    mean(R.stand.0.75[,1]), mean(R.stand.3.75[,1]), 
                    mean(R.stand.all[,1]), mean(R.short.ra[,1])),4),
  ac.m    = round(ac.times.m,2),
  s.e.m   = round(sqrt(sdv^2/(4*c(cnt,cnt,cnt,cnt,cnt,cntsc.ra)/ac.times.m)),4)
))

print(results)


# Function to compute fractions of updates that are copied/rejected, separately
# for each sort of short-cut sequence.

cr.frac <- function (cr)
{
  apply (matrix(cr,4,length(cr)/4), 1, mean)
}


# Function to plot fraction copied/rejected versus first component of state.

cr.plot <- function (R, cr, skip=0)
{
  w <- seq(1,length(cr),by=4*skip+1)
  v <- R[w,1]
  cr <- cr[w]

  plot (v,cr,type="n",xlab="",ylab="",ylim=c(0,1))

  for (i in 1:4)
  { s <- seq(i,length(cr),by=4)
    points(v[s],cr[s],pch=20, 
           col=c("darkred","darkgreen","darkblue","black")[i])
  }
}

cr.plot2 <- function (R, cr, skip=0, xlab=rep("",4))
{
  w <- seq(1,length(cr),by=4*skip+1)
  v <- R[w,1]
  cr <- cr[w]

  par(mfrow=c(2,2),mar=c(3.3,3,1,1)+0.1,mgp=c(1.8,0.7,0))

  for (i in 1:4)
  { s <- seq(i,length(cr),by=4)
    plot(v[s],cr[s],pch=20,,xlab=xlab[i],ylab="",
         ylim=c(0,1),xlim=c(-9,9),yaxs="i",xpd=NA)
  }
}


# Produce some Postscript plots.

postscript("test3stand.ps",pointsize=10,width=6.5,height=8.5,horiz=F)

par(mfrow=c(4,1),mar=c(4,3,2,1.1)+0.1,mgp=c(2,0.7,0))

w <- seq(1,nrow(R.stand.0.03),by=8)
plot(w,R.stand.0.03[w,1],
     pch=20,xaxs="i",xpd=NA,xlim=c(0,nrow(R.stand.0.03)),ylim=c(-10,10),
     ylab="",xlab="Stepsize of 0.03")
plot(w,R.stand.0.15[w,1],
     pch=20,xaxs="i",xpd=NA,xlim=c(0,nrow(R.stand.0.03)),ylim=c(-10,10),
     ylab="",xlab="Stepsize of 0.15")
plot(w,R.stand.0.75[w,1],
     pch=20,xaxs="i",xpd=NA,xlim=c(0,nrow(R.stand.0.03)),ylim=c(-10,10),
     ylab="",xlab="Stepsize of 0.75")
plot(w,R.stand.3.75[w,1],
     pch=20,xaxs="i",xpd=NA,xlim=c(0,nrow(R.stand.0.03)),ylim=c(-10,10),
     ylab="",xlab="Stepsize of 3.75")

dev.off()

postscript("test3ms.ps",pointsize=10,width=6.5,height=8.5,horiz=F)

par(mfrow=c(4,1),mar=c(4,3,2,1.1)+0.1,mgp=c(2,0.7,0))

w <- seq(1,nrow(R.stand.all),by=8)
plot(w,R.stand.all[w,1],
     pch=20,xaxs="i",xpd=NA,xlim=c(0,nrow(R.stand.all)),ylim=c(-10,10),ylab="",
     xlab="Standard Metropolis with stepsizes of 0.03, 0.15, 0.75, and 3.75")
w <- seq(1,nrow(R.short.ra),by=16)
plot(w,R.short.ra[w,1],
     pch=20,xaxs="i",xpd=NA,xlim=c(0,nrow(R.short.ra)),ylim=c(-10,10),ylab="",
     xlab="Short-cut Metropolis with stepsizes of 0.03, 0.15, 0.75, and 3.75")

dev.off()

postscript("test3c.ps",pointsize=9,width=6.5,height=4.0,horiz=F)
cr.plot2 (R.short.ra, copied.short.ra, skip=4, 
          xlab=c("Stepsize of 0.03","Stepsize of 0.015",
                 "Stepsize of 0.75","Stepsize of 3.75"))
dev.off()