Question

Hi

I need the R codes of the below Exercise
2.18 from Introducing Monte Carlo Methods with R

Answer 1

R code:

F1=function(x){
exp(-(x^2)/2) * ((sin(6*x))^2 + 3*(cos(x))^2*(sin(4*x))^2 + 1)
}
curve(F1)

$\text{approximate value of } \math{M} \text{:} \newline \text{code: integrate(F1, lower=-Inf, upper=Inf)} \newline \text{It gives : 5.89434 with absolute error less than 0.00047} \newline \text{Therefore, the normalizing constant is } \math{1/5.89434=0.1696543}$

R code:

accept_reject_alg <- function(m)
{
D <- function(x) ((sin(6*x))^2 + 3*(cos(x))^2*(sin(4*x))^2 + 1)/5
random_sample = c()
iter_num = 0
while(length(random_sample) < m){
iter_num = iter_num+1
y = rnorm(n=1,mean=0,sd=1)
u=runif(n=1,min=0,max=1)
if(u < D(y))
random_sample = c(random_sample, y)
}
RETURN = list(random_sample, iter_num)
names(RETURN) = c("sample","no_of_iterations")
return(RETURN)
}

RESULT = accept_reject_alg(2500)
RESULT$sample
RESULT$no_of_iterations

(c)

$\math{k_0 5\sqrt{2\pi } = M} \text{; and } \math{M} \text{ is the average number of iterations required to get one sample.} \newline \therefore M \approx \text{no. of iterations / sample size } \math{=5276/2500=2.1104} \newline \therefore \math{k_0\approx 0.1683856}$

code: hist(RESULT$sample)