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

k = 5
M = 10^5;
B = M/10;
h = function(x) { x }
numacc = 0;
xlist = rep(0,M);
X = 2*runif(1);  # somewhat-overdispersed starting distribution

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;
}

cat("acceptance rate =", numacc/M, "\n");
cat("mean of h is about", mean(h(xlist[(B+1):M])), "\n");

# plot(xlist[1:10000], type='l')
plot(xlist, type='l')