One of my students wrote the following code for his R exam, trying to do accept-reject simulation (of a Rayleigh distribution) and constant approximation at the same time:
fAR1=function(n){ u=runif(n) x=rexp(n) f=(C*(x)*exp(-2*x^2/3)) g=dexp(n,1) test=(u<f/(3*g)) y=x[test] p=length(y)/n #acceptance probability M=1/p C=M/3 hist(y,20,freq=FALSE) return(x) }
which I find remarkable if alas doomed to fail! I wonder if there exists a (real as opposed to fantasy) computer language where you could introduce constants C and only define them later… (What’s rather sad is that I keep insisting on the fact that accept-reject does not need the constant C to operate. And that I found the same mistake in several of the students’ code. There is a further mistake in the above code when defining g. I also wonder where the 3 came from…)