Monte Carlo integration method: sampling in region?

[Curl]

I want to use this method to evaluate a triple integral that is otherwise impossible (not elementary). The region is the unit sphere.

The Monte Carlo method sounds simple, but if I were to use MATLAB to do this, how would I pick my random sample points so that they are INSIDE the region?

Are there any efficient methods for doing this, besides chosing points randomly from the unit cube and checking them all to see that they are inside the sphere? This selection process sounds inefficient, especially for regions other than spheres.

Does anyone have any tips?

 

[Office_Shredder]

You could just describe the point in polar coordinates and pick random points in the space:
radius lies inside [0,1]
theta in [0,2pi]
phi in [0,pi]

In general your ability to pick random points is only as good as your ability to describe the region to the computer

 

[CRGreathouse]

You could just describe the point in polar coordinates and pick random points in the space:
radius lies inside [0,1]
theta in [0,2pi]
phi in [0,pi]

If you pick the radius uniformly at random from [0, 1] you pick near the center too much.

 

[jackmell]

I want to use this method to evaluate a triple integral that is otherwise impossible (not elementary). The region is the unit sphere.

The Monte Carlo method sounds simple, but if I were to use MATLAB to do this, how would I pick my random sample points so that they are INSIDE the region?

Are there any efficient methods for doing this, besides chosing points randomly from the unit cube and checking them all to see that they are inside the sphere? This selection process sounds inefficient, especially for regions other than spheres.

Does anyone have any tips?

Would you not pick the points randomly from the range (-1,1)? That is, from a cube 2x2x2? centered at the origin? Then check if they're in the unit sphere, then average the value of the function over the points inside the sphere and then multiply by the volume of the sphere? Not real familiar with the Monte Carlo method but I'm a programmer and thought I'd try that in Mathematica (don't have Matlab). The code below integrates the function f(r)=r^2-r^3 over the unit sphere via Monte Carlo method with 10000 points and then computes the integral by more standard means. You probably already know how to code it. Just thought I'd do it for fun. The results look pretty close to the explicit numerical result.

In[70]:=
myf[r_] := r^2 - r^3; 
mymax = 10000; 
mysum = 0; 
mynet = 0; 
maxr = 1; 
For[i = 1, i <= mymax, i++, 
   mypt = {RandomReal[{-maxr, maxr}], 
      RandomReal[{-maxr, maxr}], 
      RandomReal[{-maxr, maxr}]}; 
    mynorm = Norm[mypt]; 
    If[mynorm < maxr, mynet++; 
      mysum += myf[mynorm]; ]; ]; 
myavg = mysum/mynet; 
myint = (4/3)*Pi*maxr^3*myavg
N[8*Integrate[r^2*myf[r]*Sin[\[Phi]], 
    {\[Phi], 0, Pi/2}, {\[Theta], 0, Pi/2}, 
    {r, 0, maxr}]]

Out[77]=
0.4185168995775403

Out[78]=
0.41887902047863906

 

[CRGreathouse]

Are there any efficient methods for doing this, besides chosing points randomly from the unit cube and checking them all to see that they are inside the sphere? This selection process sounds inefficient, especially for regions other than spheres.

For a 3-sphere, it's not too bad. For higher dimensions, it gets progressively worse. By dimension 50, it's essentially impossible to use this method (the 50-sphere's content is something like d^50 / 6e27, where d is the diameter).

 

[mathman]

You can use spherical coordinates, but make sure you use the coordinates properly.
dV=r²drdθsinφdφ. Therefore you need to pick r³ uniform - domain (0,1), θ uniform - domain (0,2π), and cosφ uniform - domain (-1,1).

 

[Curl]

You can use spherical coordinates, but make sure you use the coordinates properly.
dV=r²drdθsinφdφ. Therefore you need to pick r³ uniform - domain (0,1), θ uniform - domain (0,2π), and cosφ uniform - domain (-1,1).

If I use that change of variables the integrand becomes extremely messed-up.

If you are suggesting this for Monte Carlo, then like someone else said, there will be too many points chosen towards the core of the sphere (if you randomly select values for r) which will screw up the answer.

I guess I'll try to select random points from -1 to 1 for x, then for y I would use the constraint of the unit circle in the xy plane, then for z use the third constraint. I don't know how to use mathematica so I'll try this in MATLAB, though your answer looks good jackmell.

 

[Office_Shredder]

If you are suggesting this for Monte Carlo, then like someone else said, there will be too many points chosen towards the core of the sphere (if you randomly select values for r) which will screw up the answer.

He's suggesting that you use that to get the probability distribution from which you pick the radius and the angles

 

[mathman]

If I use that change of variables the integrand becomes extremely messed-up.

If you are suggesting this for Monte Carlo, then like someone else said, there will be too many points chosen towards the core of the sphere (if you randomly select values for r) which will screw up the answer.

I guess I'll try to select random points from -1 to 1 for x, then for y I would use the constraint of the unit circle in the xy plane, then for z use the third constraint. I don't know how to use mathematica so I'll try this in MATLAB, though your answer looks good jackmell.

If you pick from the distributions for r and the angles as I suggested, you will get points uniformly distributed in the volume of the unit ball. For calculation purposes it is very easy to get x, y, and z, which then can be used for sampling the integrand.

 

[Dickfore]

Seeing your region of integration is a unit sphere, I suspect you can go over to spherical coordinates. If you have a uniform distribution inside a sphere of radius R, i.e. if the distribution is:

[tex]
\varphi(x, y, z) = \frac{3}{4 \, \pi \, R^{3}} \, \Theta(R - \sqrt{x^{2} + y^{2} + z^{2}})
[/tex]

where

[tex]
\Theta(x) = \left\{\begin{array}{l}{1, \; x \ge 0 \\ 0, \; x < 0 \end{array}\right.
[/tex]

is the Heaviside step function and you know the Jacobian to be:

[tex]
J = \frac{\partial(x, y, z)}{\partial(r, \theta, \phi)} = r^{2} \, \sin{\theta} \, dr \, d\theta \, d\phi
[/tex]

then:

[tex]
\varphi(x,y,z) \, dx \, dy \, dz = \tilde{\varphi}(r, \theta, \phi) dr \, d\theta \, d\phi
[/tex]

and you express everything in spherical coordinates, you will get:

[tex]
\tilde{\varphi}(r, \theta, \phi) = \frac{3}{R^{3}} \theta(R - r) \, r^{2} \, \frac{1}{2} \, \sin{\theta} \, \frac{1}{2 \pi}
[/tex]

You can see that this distribution function factors out into products of function which depend on a single variable. This means you can obtain the same distribution (uniform inside a sphere of radius R) if you generate two one-dimensional random variables with the following distributions:

1. [tex]\varphi_{1}(\phi) = \frac{1}{2 \, \pi}, \; 0 \le \phi \le 2 \pi[/tex];

2. [tex]\varphi_{2}(\theta) = \frac{1}{2} \, \sin{\theta}, \; 0 \le \theta \pi[/tex];

3. [tex]\varphi_{3}(r) = \frac{3 r^{2}}{R^{3}}, \; 0 \le r \le R[/tex]

and, once you have a particular sample [itex](r_{0}, \theta_{0}, \phi_{0})[/tex], then you can calculate the Cartesian coordinates:

[tex]
x_{0} = r_{0} \, \sin{\theta_{0}} \, \cos{\phi_{0}}
[/tex]
[tex]
y_{0} = r_{0} \, \sin{\theta_{0}} \, \sin{\phi_{0}}
[/tex]
[tex]
z_{0} = r_{0} \, \cos{\theta_{0}}
[/tex]

It is guaranteed that the sample of points [itex](x_{0}, y_{0}, z_{0})[/itex] is uniformly distributed inside the sphere of radius R.

The advantage that this method has is that there are no points that are no sample points being rejected. The ratio of the volumes of a sphere with radius R and a cube with side 2R (so that the sphere lies wholly inside the cube) is:

[tex]
\frac{\frac{4 \, \pi \, R^{3}}{3}}{(2 \, R)^{3}} = \frac{\pi}{6} = 0.52
[/tex]

so, by simply rejecting points, you will need to generate approximately twice as many as you need.

The disadvantage is that you need generate random variables with non-uniform distribution. Fortunately, the distributions are elementary and you can use the following Theorem:

If [itex]f(x), \alpha \le x \le \beta[/itex] is a continuous function that represents the probability distribution function of a random variable X then:

[tex]
dY = f(X) dX \Rightarrow Y = F(X) = \int_{\alpha}^{x}{f(t) \, dt}, \alpha \le x \le \beta
[/tex]

the random variable Y has uniform distribution on the segment [itex][0, 1][/itex]. If you can invert this relationship:

[tex]
X = F^{-1}(Y), 0 \le y \le 1
[/tex]

and you generate a uniformly distributed random variable, then X will have the PDF required. In our case:

[tex]
Y = \int_{0}^{r}{\frac{3 x^{2}}{R^{3}} \, dx} = \left(\frac{r}{R}\right)^{3} \Rightarrow r = R \, Y^{1/3}
[/tex]

[tex]
Z = \int_{0}^{\theta}{\frac{\sin{t} \, dt}{2}} = -\left.\frac{\cos{t}}{2}\right|^{\theta}_{0} = \frac{1 - \cos{\theta}}{2} = \sin^{2}{\frac{\theta}{2}} \Rightarrow \theta = 2 \, \arcsin{Z^{1/2}}
[/tex]

These formula furnish the whole algorithm for generating points uniformly inside a sphere of radius R. Notice, however, that Monte Carlo is much more efficient if the points are generated with a distribution that is comparable to the values of the integrand.

 

[Dickfore]

I have generated 10,000 points in less than 5 seconds using the above algorithm in Mathematica.

 

[Curl]

Very nice. For a sphere I would suspect that a distribution equation be elementary. However if I want to generalize this method for any integration limits, I'd have to do something different.

I might try to make a MATLAB function in case I need to use integral approximations extensively. I think the rejection method is the simplest to program, even if lots of points are rejected.

Isn't it still quicker (computation time not efficiency) to generate random numbers and reject those that do not match? Or at least do a procedure like this, example for a sphere:

Chose a random number for X between -1 and 1.
For each X, chose a random number for Y such that y < sqrt(1-x^2)
For each X and Y, chose a random number for Z such that Z < sqrt(1-X^2-Y^2)
Repeat this every time you want a set of (x,y,z) coordinates (points).

I don't know if this is any quicker or easier to generalize than simply choosing points in a box containing the region and rejecting them afterward.

 

[Dickfore]

Do you know what the Monte Carlo Method's strength actually is?

 

[Curl]

For me the strength is that it allows me to evaluate integrals that are not elementary or that would take me too long to work out by hand.

 

[Dickfore]

There are much more efficient methods for numerical integration, though. Why Monte Carlo?

 

[mathman]

There are much more efficient methods for numerical integration, though. Why Monte Carlo?

This is true for one or two dimensional integrals. However Monte Carlo methods are better for integrals with more than two dimensions.