[Math] Estimate $\pi$ using the Monte Carlo Method in MATLAB

MATLABmonte carlonumerical methods

Evaluate the area of a circle of radius $1= \pi$ using Monte Carlo method . Hence we can generate pairs of random numbers $(x_i,y_i) \in [-1,1]$.

Thus :

$$ \pi= \frac {Number Of Samples Inside The Circle}{Total Number Of Samples} X 4 $$

Consider a unit circle inscribed in a square, each of the small circles drawn on this figure represents a random point that was generated in the square, the red and blue circles represent points inside and outside the unit circle respectively (I got the figure from google) . If we choose a point uniformaly at random within the square, then the probability that the point lies in the unit circle is

$$ \pi= \frac {Number Of Samples Inside The Circle}{Total Number Of Samples} $$

We know that the area of the circumscribed square is $4$, if we knew $p$, then we could compute the area of the unit circle:

Area of the unit circle = $p$ x area of the circumscribed square = $4p$

How to Estimate $\pi$ using the Monte Carlo Method in MATLAB, then what's the error a quarter of a circle ?

Best Answer

The following code is vectorized, and thus typically faster than using loops:

N = 1e6;                % Number of samples
r = 1;                  % Circle radius
x = 2*rand(1,N)-1;      % N samples between -1 and 1
y = 2*rand(1,N)-1;      % N samples between -1 and 1
r2 = x.^2 + y.^2;       % compute squared distance to origin
area = mean(r2<=r^2)*4; % the area is 4 (area of square) times the proportion
                        % of points in the circle. Note that the threshold
                        % should be r^2 (1^2)
error = pi - area       % error in the estimation of pi

N_all = round(10.^(1:.5:7));      % Numbers of samples
r = 1;                            % Circle radius
errors = NaN(size(N_all));        % Preallocate result
for k = 1:numel(N_all)            % Loop over size of N_all
    N = N_all(k);                 % Current N
    x = 2*rand(1,N)-1;            % N samples between -1 and 1
    y = 2*rand(1,N)-1;            % N samples between -1 and 1
    r2 = x.^2 + y.^2;             % compute squared distance to origin
    area = mean(r2<=r^2)*4;       % the area is 4 (area of square) times the
                                  % proportion of points in the circle
    error_all(k) = pi-area;       % error in the estimation of pi
end                               % End loop
semilogx(N_all, error_all, 'o-'); % Plot with logarithmic x axis
grid                              % Use grid

Related Solutions

[Math] Using Monte Carlo method to estimate the value of $\int_0^1\int_0^1e^{xy} \ dx \ dy$

Here is a code that should do the trick:

s=0;
for i=1:N  
x=rand(0,1); 
y=rand(0,1);
s=s+exp(xy); 
end 
I = s/N;

[Math] Monte Carlo gamma function

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.

Pick $n$ randomly distributed points $x_1, x_2, \ldots, x_n$ in the interval $a \le x \le b$.
Determine the average value of the function $\displaystyle f_{avg} = \dfrac{1}{n} \sum_{i = 1}^n f(x_i)$.
Compute the approximation to the integral $\displaystyle \int_a^b f(x)~dx = (b-a)f _{avg}$.
An estimate of the error is Error $\approx (b-a) \sqrt{\dfrac{f_{avg}^2 - (f_{avg})^2}{n}}$, where $f_{avg}^2 = \dfrac{1}{n} \sum_{i = 1}^n f^2(x_i)$.

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.

First, we must determine the rectangular $R$ box containing the area $A$ where $R = \{ (x, y): a \le x \le b ~\mbox{and}~ 0 \le y \le d$ where $\displaystyle d = \max_{a \le x \le b} f(x)$.
Second, randomly pick points $n$ point, $(x_n, y_n) \in R$, where the $x_n$ and $y_n$ are chosen from an independent uniformly distributed random variables over $x_n \in [a, b]$ and $y_n \in [0, d]$ , respectively.
Third, calculate the ratio $\rho$ as $\rho = \dfrac {m}{n}$, where $m$ is number of points that lie in $A$.
The area is computed using the approximation area $A = \rho (b-a) d$.
An estimate for the accuracy of the above computation is Error $\approx (b-a) d \sqrt{\dfrac{\rho - \rho^2}{n}}$.

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:

$~~~~~~~~~~~~~$ enter image description here

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:

$~~~~~~~~~~~~~$ enter image description here

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$$

Best Answer

Related Solutions

[Math] Using Monte Carlo method to estimate the value of $\int_0^1\int_0^1e^{xy} \ dx \ dy$

[Math] Monte Carlo gamma function

Related Question