## code to accompany the importance samplin example from
## lecture 6 of PSTAT 215A, Bayesian Inference, UCSB
##
## Robert B. Gramacy

## first example, based on the Cauchy distribution

## rh:
##
## draw from h=2/x^2 by the method of inversion

rh <- function(n=1)
  {
    u <- runif(n)
    return(2/(1-u))
  }

## plot a cartoon of the IS setup
xt <- seq(-6,6,length=1000)
plot(xt, dt(xt, df=1), type="l", lwd=2, ylim=c(0,0.45),
     xlab="theta", ylab="f(theta)", bty="n")
abline(v=2, col=2, lty=2, lwd=2)
xt2 <- xt[xt >= 2]
lines(xt2, 2/(xt2^2), lty=3, col=3, lwd=3)
legend("left", c("Cauchy", "theta=2", "h(theta)"),
       lty=1:3, col=1:3, lwd=2, bty="n")

## take n samples from the g distribution
n <- 1000
y <- rh(n)

## calculate the importance weights
w <- y^2/(2*pi*(1+y^2))

## calculate the IS estimator est-is
est.is <- sum(w)/n

## calculate the plug-in estimator based on the
## same number of samples
x <- rt(n, df=1)
est.plugin <- sum(x > 2)/length(x)

## calculate the truth
truth <-  pt(2, df=1, lower.tail=FALSE)

## compare the two estimates to the truth
est <- data.frame(hat=est.is, plugin=est.plugin,
                    truth=truth)
est

## don't need to go back to slides here

## look at the variance of the IS and plugin estimators
T <- 1000
hat <- plugin <- rep(NA, T)
for(t in 1:T) {
  y <- rh(n)
  w <- y^2/(2*pi*(1+y^2))
  hat[t] <- sum(w)/n
  x <- rt(n, df=1)
  plugin[t] <- sum(x > 2)/length(x)
}
dplugin <- density(plugin)
dhat <- density(hat)
plot(dplugin, lwd=2, ylim=range(dhat$y), bty="n", main="",
     xlab="MC estimates")
lines(dhat, col=2, lty=2, lwd=2)
abline(v=truth, col=3, lty=3, lwd=2)
legend("right", c("plug-in", "IS", "truth"), col=1:3,
       lty=1:3, lwd=2, bty="n")


## second example, illustrating an estimate of the
## CDF and sample importance resampling for the Beta
## distribution

## take n proposals from the g (uniform) distribution
n <- 1000
y <- runif(n)

## calculate the importance weights and then normalize
## them
w <- dbeta(y, 3, 2)/dunif(y)
w <- w/sum(w)

## create the resulting emperical CDF
o <- order(y)
cdf <- rep(0, n)
cdf[1] <- w[o[1]]
for(i in 2:n) cdf[i] <- cdf[i-1] + w[o[i]]

## plot the emperical cdf versus the truth
xt <- seq(0,1,length=1000)
plot(xt, pbeta(xt, 3, 2), type="l", lwd=3, xlab="x",
     ylab="F(x)", lty=2, bty="n")
lines(y[o], cdf, type="l", lwd=3, col="blue")

## take a SIR sample according to the weights
x <- sample(y, n/2, replace=TRUE, prob=w)
hist(x, lwd=2, lty=2, freq=FALSE, breaks=20, bty="n", main="")
lines(xt, dbeta(xt, 3, 2), lwd=2)