Sample Code/Rejection_Sampling_Example.R

In jstarling1/starlib: All of my useful functions.

#MCMC Exercise 2: Rejection Sampling

#Sample from f(x) proportional to: x^3 * exp(-2x) * (1-exp(-x))^2.

#Let h(x) = x^3 * exp(-2x), which has form of gamma(shape=4, rate=2).
#Let l(x) = (1-exp(-x))^2.

#Sample x~h(x) and y~U(0,1) independently.


#Let M = 3/8.
M = 3/8

#Set up lx function for convenience.
lx_fun = function(x){
	return((3/8)*(1-exp(-x))^2)
}

#------------------------------------------------------

#Repeat B times to generate a good sample of fx,
#and to approximate the acceptance probability.

B = 10000		#Number of iterations.
fx = NA			#Empty vector to hold fx realizations.
accepts = 0		#Hold number of accepts.

for (i in 1:1000){	
	M = 3/8				#Upper bound for function lx, found by taking limit as x->inf.
	y = runif(1,0,1)	#Generate a U(0,1) 
	hx = rgamma(1,shape=4,rate=2)	#Sample from gamma dist.
	lx = lx_fun(hx)		#Plug hx into lx.

	bound = lx / M		#Calculate bound for accept/reject.

	#If bound greater than the uniform realization, accept hx as a realization of fx:
	if(y <= bound) {
		fx = c(fx,hx)
		accepts = accepts + 1
	}
}

#Probability of acceptance.
p_accept = sum(decisions)/B
p_accept

#hx is the sample from fx.

#------------------------------------------------------
#PLOTTING:

fx_fun = function(x){
	return(x^3 * exp(-2*x) * (1-exp(-x))^2)
}

hist(fx,probability=T)
lines(x=seq(0,10,by=.01),y=4*fx_fun(seq(0,10,by=.01)),col='blue')