Question

Q. 2 a) Using only Unif 0,1) random variates, use a Monte Carlo algorithm to approximate the value of the Gamma function 0 by considering the function as an expectation of a function of a random variable. b) Show that if a Monte Carlo simulation of size N is used then the variance of the Monte Carlo estimator is Var (「(a))- 2a-1)-[「(a)]2] provided that α > 0.5. c) Write an R function to implement the method returning the estimated value of the function and the standard error of the estimation. Examine how well the method works for various values of α and various simulation sizes.

Answer 1

0 Then ri 2 a -,

R Program

f=function(n,alpha)
{
s1=0
s2=0
x=numeric(n)
for(i in 1:n)
{
x[i]=-log(1-runif(1,0,1))
s1=s1+x[i]^(alpha-1)
s2=s2+x[i]^(2*alpha-2)
}
est=s1/n
var_est=s2/n-est^2
se_est=sqrt(var_est/n)
cat(" ",c("Estimated=",est,"Actual=",gamma(alpha),"SE=",se_est))
}
# Values for Gamma(3)
f(100,3)
f(1000,3)
f(5000,3)
f(10000,3)

#Values for Gamma(6)
f(100,6)
f(1000,6)
f(5000,6)
f(10000,6)

R Output

> f(100,3)

Estimated= 2.39643722309558 Actual= 2 SE= 0.78419283402655> f(1000,3)

Estimated= 2.39853899997827 Actual= 2 SE= 0.184532664126924> f(5000,3)

Estimated= 1.91312239062837 Actual= 2 SE= 0.059122563425076> f(10000,3)

Estimated= 2.04307398594427 Actual= 2 SE= 0.0475782103374921>
> #Values for Gamma(6)
> f(100,6)

Estimated= 69.9108955051665 Actual= 120 SE= 50.9203707580712> f(1000,6)

Estimated= 109.469711149383 Actual= 120 SE= 23.651514011367> f(5000,6)

Estimated= 154.31873737748 Actual= 120 SE= 33.2707135100746> f(10000,6)

Estimated= 120.553249455824 Actual= 120 SE= 17.2548285523348

Comment

For fixed alpha, as we increase simulation number, the monte carlo estimate becomes closer to the actual with a substantial decrease in SE value. Only alpha=3 and 6 are considered here but works well for other choices also.