# independence sampler example with pi = Exp(1), q = Exp(k)

k = 5
M = 10^4;
B = M/10;
h = function(x) { x }
numacc = 0;
xlist = rep(0,M);
X = rexp(1);  # starting distribution equals pi

for (i in 1:M) {
    Y = rexp(1,k);
    A = exp((k-1)*(Y-X));
    U = runif(1);
    if (U<A) {
	# accept proposal
	X = Y;
	numacc = numacc + 1;
    }
    xlist[i] = X;
}

hlist = xlist  # Here "h" is just the identity functional.

cat("Independence sampler: run length", M, "with burn-in", B, "\n");
cat("acceptance rate =", numacc/M, "\n");
estmean = mean(hlist[(B+1):M])
cat("mean of h is about", estmean, "\n")

se1 =  sd(hlist[(B+1):M]) / sqrt(M-B)
cat("iid standard error would be about", se1, "\n")

varfact <- function(xxx) { 2 * sum(acf(xxx, plot=FALSE)$acf) - 1 }
se2 = se1 * sqrt( varfact(hlist) )
cat("true standard error is about", se2, "\n")

cat("95% confidence interval is about (",
      estmean - 1.96*se2, ",", estmean + 1.96*se2, ")\n")

# plot(xlist, type='l')