We are trying to use Monte Carlo Simulation to find:
$$\displaystyle \int_0^\infty \dfrac 34 x^4e^{-x^{3/4}}\,dx$$
As you have discovered, the infinite limits are problematic. We can do several things to overcome this
Approach 1: A naive approach is to plot the function given that it has a decaying exponential which tends to zero very quickly. If we plot this function, we see:
So, lets try to approximate it as:
$$\displaystyle \int_0^{50} \dfrac 34 x^4e^{-x^{3/4}}\,dx$$
Using your approach (recall, now your code needs to change to generate random numbers $(a, b)$ and with the function of interest) with $n = 10000$, we get:
$$I = 387.256$$
If we calculate the exact integral we get:
$$I = \Gamma \left(\frac{20}{3}\right) \approx 389.035$$
Approach 2: In this approach we will make a change in variables so as to map an infinite range of integration into a finite one.
A good example that works for functions that decrease towards infinity faster than $\frac{1}{x^2}$ is to use:
$$\displaystyle \int_a^b f(x)~dx = \int_{\frac{1}{a}}^{\frac{1}{b}} \dfrac{1}{t^2}f\left(\dfrac{1}{t}\right)~dt$$
With this altered approach, we can now do integrals when one of the limits is infinite
and the other limit is of the same sign. If you need to integrate say something that has $a = 0$, what you do is break the integral into two parts and use Monte Carlo on each one.
If we have infinite limits for both limits, we could split the integral into three or more parts. So, for this problem, we have to split it into two pieces as:
$$\displaystyle \int_0^\infty \dfrac 34 x^4e^{-x^{3/4}}\,dx = \int_0^1 \dfrac 34 x^4e^{-x^{3/4}}\,dx + \int_0^1 \dfrac{1}{t^2} \left(\dfrac 34\right) \left(\dfrac 1t \right)^4 e^{{-(\frac 1t)}^{\frac 34}}~dt$$
Using the Monte Carlo for the first part, we get $I_1 = 0.0647226$ and for the second part, we get $I_2 = 387.551$, which yields:
$$ \displaystyle \int_0^\infty \dfrac 34 x^4e^{-x^{3/4}}\,dx = 387.616$$
Compare that to our previous and the exact result.
Approach 3: In this approach we will make a change in variables so as to map an infinite range of integration into a finite one. This time, we will transform variables using $u = \frac{1}{1+x}$, yielding:
$$\displaystyle \int_0^\infty \dfrac 34 x^4e^{-x^{3/4}}\,dx = \int_0^1 \frac{3 e^{-\left(\frac{1-u}{u}\right)^{3/4}} (1-u)^4}{4 u^6} ~ du$$
Running the Monte Carlo simulation with $n = 10000$ yields:
$$I = 390.18704821709946$$
Recall that if $Y$ is a random variable with density $g_Y$ and $h$ is a bounded measurable function, then $$\mathbb E[h(Y)] = \int_{\mathbb R} h(y)g_Y(y)\,\mathsf dy. $$
Moreover, if $Y\sim\mathcal U(0,1)$, then $a+(b-a)U\sim\mathcal U(a,b)$. So applying the change of variables $x=a+(b-a)u$ (with $a=0$, $b=\pi$) to the given integral, we have
$$I = \int_0^1 \frac{\pi}{\sqrt{2\pi}} e^{-\frac12\sin^2 (\pi u) }\,\mathsf du=\int_0^1 h(u)\,\mathsf du, $$
with $h(u)=\sqrt{\frac\pi 2} e^{-\frac12\sin^2 (\pi u) }$. It follows then that $I=\mathbb E[h(U)]$ with $U\sim\mathcal U(0,1)$. Let $U_i$ be i.i.d. $\mathcal U(0,1)$ random variables and set $X_i=h(U_i)$, then for each positive integer $n$ we have the point estimate
$$\newcommand{\overbar}[1]{\mkern 1.75mu\overline{\mkern-1.75mu#1\mkern-1.75mu}\mkern 1.75mu}
\widehat{I_n} =: \overbar X_n= \frac1n \sum_{i=1}^n X_i$$ and the approximate $1-\alpha$ confidence interval
$$\overbar X_n\pm t_{n-1,\alpha/2}\frac{S_n}{\sqrt n}, $$
where $$S_n = \sqrt{\frac1{n-1}\sum_{i=1}^n \left(X_i-\overbar X_n\right)^2} $$ is the sample standard deviation.
Here is some $\texttt R$ code to estimate an integral using the Monte Carlo method:
# Define "h" function
hh <-function(u) {
return(sqrt(0.5*pi) * exp(-0.5 * sin(pi*u)^2))
}
n <- 1000 # Number of trials
alpha <- 0.05 # Confidence level
U <- runif(n) # Generate U(0,1) variates
X <- hh(U) # Compute X_i's
Xbar <- mean(X) # Compute sample mean
Sn <- sqrt(1/(n-1) * sum((X-Xbar)^2)) # Compute sample stdev
CI <- (Xbar + (c(-1,1) * (qt(1-(0.5*alpha), n-1) * Sn/sqrt(n)))) # CI bounds
# Print results
cat(sprintf("Point estimate: %f\n", Xbar))
cat(sprintf("Confidence interval: (%f, %f)\n", CI[1], CI[2]))
For reference, the value of the integral (as computed by Mathematica) is $$e^{-\frac14}\sqrt{d\frac{\pi }{2}} I_0\left(\frac{1}{4}\right) \approx 0.991393, $$
where $I_\cdot(\cdot)$ denotes the modified Bessel function of the first kind, i.e. $$I_0\left(\frac14\right) = \frac1\pi\int_0^\pi e^{\frac14\cos\theta}\,\mathsf d\theta. $$
Best Answer
Fixing $a, b, c$, you should be evaluating points at
u_xandu_yrather thanxandy.Formula for
u_xandu_yareunif(10000).Edit: