This question was asked before but I'd like to ask something more precise given the answer that was given.
[ Estimate gamma function using monte carlo ]
What is the criterion for a random point to be considered inside the function gamma?
Usually the Monte Carlo method multiplies the area of a square enclosing the function by the number of points generated ending up inside the function, divided by the total number of points generated.
Here, I'm confused as to what is the area of the square enclosing gamma and what is the number of points generated ending up inside the function. I'm guessing (1/N) is the division by the total number of points generated. To me, it seems like this is just an estimate of the gamma function itself using a sum, but not actually a Monte Carlo method, so I might be missing something.
Anyone could help? Thank you !
Best Answer
This is a very good question.
For one-dimensional integrals, there are two ways we could approach it.
Method 1: This is the approach in the other link you listed.
Notice that in this approach, we know the function values apply because we are using the function itself as part of the calculation. Just like other numerical integration techniques, you can interpret $(1/n)f(x_n)$ as the area of a rectangle of height $f(x_n)$ and base $1/n$, so you are adding up the areas of those rectangles. Typically, the more rectangles (larger the $n$), the better the approximation.
Method 2: The way you are asking.
It is worth noting that Method 2 makes it easier to move to higher dimensional integrals, but is probably overkill for $1D$. Lets apply these methods and compare to the exact problem you mentioned for the Gamma Function in the other problem:
$$\displaystyle \Gamma(\beta) = \int_0^\infty x^{\beta - 1} e^{-x} dx$$
I will use $\beta = \pi$.
Approach 1: A naive approach is that since we have this integral going to infinity, we can look at a plot and determine a suitable range to do the Monte Carlo simulation.
$$\displaystyle \Gamma(\pi) = \int_0^\infty x^{\pi - 1} e^{-x} dx$$
If we plot this function over $0 \le x \le 10$, we get a reasonable range given how fast the exponential approaches zero:
$~~~~~~~~~~~~~$
Using Method 1, with $n = 1000$, we get:
$$I = 2.28421$$
Using Method 2, with $n = 1000$, we get:
$$I = 2.29143$$
Calculating the exact value, we have:
$$\Gamma[\pi] = 2.288037795340033$$
Both approximations are not bad for only $1000$ samples and more samples will improve the result given the Central Limit Theorem. Note running multiple trials does show that there is some variation in the result, but in all cases it is not bad.
Lastly, using a sophisticated Adaptive-Quasi-Monte-Carlo Method, the samples and plot look like:
$~~~~~~~~~~~~~$
This approach yields a result of:
$$\Gamma[\pi] = 2.28804242000413$$
Update
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 \Gamma(\pi) = \int_0^\infty x^{\pi - 1} e^{-x} dx = \int_0^1 x^{\pi - 1} e^{-x} dx + \int_0^1 \dfrac{1}{t^2} \left(\dfrac 1t\right)^{\pi - 1} e^{-\frac 1t} ~dt$$
Using the Monte Carlo for the first part, we get $I_1 = 0.14972956088031433$ and for the second part, we get $I_2 = 2.1272262548687735$, which yields:
$$ \displaystyle \Gamma(\pi) = \int_0^\infty x^{\pi - 1} e^{-x} dx = 2.276955815749088$$
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:
$$\int_0^1 \left( \frac{1-u}{u}\right)^{\beta-1}e^{-\frac{1-u}{u}}\frac{1}{u^2}du$$
Using Method I with $n=1000$ yields:
$$I = 2.2895585648008825$$