# Monte Carlo Integration Interval Probability

Use MC integration to estimate the probability that X * exp(X) < 2.5, assuming that X ~ Gamma(1.2,3.7)

data mcprobdata;
/* generate the sample */
call streaminit(23891);
count = 0;
do i=1 to 2000;
p = rand('uniform');
x = quantile('Gamma', p, 1.2, 3.7);
count = count + (x*exp(x) < 2.5);
end;

/* calculate the estimate */
t = count/2000;
se = sqrt(t*(1-t)/2000);

run;


Is my t variable correctly calculating this interval probability?

• This probability can be found using other techniques of integration: why don't you compare those results with yours? It might help to know that $x\exp(x)\lt 2.5$ is equivalent to $x\lt 0.958586$ (approximately). – whuber Apr 18 '14 at 14:13

Solution in R for reference.

N <- 10^6
x <- rgamma(N, shape =  1.2, scale = 3.7)
sum(x * exp(x) < 2.5) / N
 0.155834


Check the SAS documentation to see how they parameterize the gamma distribution (shape/scale or shape/rate?).

Yes, it is, as long as x is successfully being drawn as intended.

2000 is quite a low number of MC draws to use.

You can test it by doing a numerical integration in Excel or checking the Gamma cdf by some other means.