High dimension Monte Carlo Integration

[Ellen Zhao]

I am implementing a Monte-Carlo-Integration in java:
Integrand: Math.exp( Math.sin(x1)*Math.sin(x2)*Math.sin(x3)*Math.sin(x4)*Math.sin(x5)*Math.sin(x6)*Math.sin(x7)*Math.sin(x8)*Math.sin(x9)*Math.sin(x10) )dx1dx2...dx10
Interval: [0, Math.PI/2]^10
Number of sample dots: 10000000
I know how to compute integration with trapeziodal algorithm and Simpson�s algorithm, and I know how to compute PI with Monte-Carlo algorithm, but how can I deal with this integrand?
I have to tell my program the acceptance condition, but how to specify the acceptance area here? Thank you very much in advance!

Regards,
Ellen
[ January 16, 2003: Message edited by: Ellen Fu ]

[Norm Miller]

Hi Ellen,
Why don't you simplify the problem by forgetting about the x2 through x10, reducing it to one dimension (and one integral)? You don't have to code it.
Then bring in the x2 variable and think that through. By that point you should have an idea of how it will work in ten dimensions.
The good news is that since sin(x) is always between 0 and 1 ( for the intervals in question ), you know that exp(sin(x)) is between 1 and e. So you know the dimensions of the bounding box.
For those who haven't followed this argument, the following page should help you understand:
[URL=http://www.taygeta.com/rwalks/node5.html[/URL]

[Ellen Zhao]

Hi Norm,
Thank you very much for your enlightenment. I�ll go back to my freshman caculus to see how to reduce the dimension for this integrand.
Regards,
Ellen

[Norm Miller]

Hi again,
All I was thinking for the first step was to think about the function e(sin(x1)) integrated over 0 to pi/2. To do that by monte-carlo technique, you just generate a lot of x's and y's randomly ( the x's constrained to 0, pi/2, the y's constrained to 1, e). Each x has a corresponding y. Then you just see (by evaluating the function) whether the y value is less than the value of the function for the corresponding x value. The integral is approximated by the ratio of points which are in the region to total points.
In one dimension, that is fairly intuitive. Now bring in x2. The technique is not too difficult to extend and it isn't hard to visualize what's happening.
When you get to x3, x4, etc, visualization (for me at least) is about impossible, but by now I think I can just proceed as before and do the math. Of course, by now I probably know how to write the program, which is the idea.
Another thing. I think you can integrate the one-dimensional case as you did in calculus class and check the answer you get if you program the one dimensional case for practice. That way you can make sure you didn't forget any scale factors or whatever. Then you will be more certain you have the multi-dimension integration program correct.
Good Luck
Norm

[Ellen Zhao]

Hi Norm,
Thank you so much! I�ll put the output here when I have completely implemented your idea.
Best Regards,
Ellen

[Ellen Zhao]

got some very buggy result-tested with 10^n( n from 2 to 7 )random numbers, modifying the code...
n = 100 I = 91.45565289510405
n = 1000 I = 91.45341951653647
n = 10000 I = 91.4531961786797
n = 100000 I = 91.45317384489402
n = 1000000 I = 91.45317161151546
n = 10000000 I = 91.45317138817761

Regards,
Ellen

[Norm Miller]

Looks like you learned that about a thousand iterations gets you as close as you need to be.
But is the result approximately correct? I found the answer (numerically) in 1 dimension to be about 3.1 I don't think the added dimensionality will change that a lot, except to make it smaller (??)
Arm-waving argument follows:
exp(sin(x1)) between 0 and pi/2 ranges from 1 to e (because sin(x1) is between 0 and 1.
So the 3.1 result is at least reasonable.
Now try exp(sin(x1)sin(x2)). Unless x2 happens to be 1, the product of the two sins is smaller than sin(x1). So I would think the integral would get smaller as the dimensionality increases.
I guess I am just suggesting that you work up an estimate of the result you should get.

[Ellen Zhao]

Hi Mr. Miller,
Yes, your approximate result is absolutely reasonable. And I understand your idea.
My classmates and I discussed this question. Actually there exists some universal formular for computing high-dimensional MC Integration. However, the text we got was quite obscure, that led a magnificent variety of our results, no one was close to the very reasonable 3.1.
Here is the algorithm:
1. n: the number of random number
k: the dimension
it asserts there would be n*k calls to the random number generator: c1,c2,...ck, ck+1,...c2k,...,cnk.
2. Compute the Yn which is approximate to I(f):
b = PI/2
a = 0
A = ( \Pi( capitalized 16th letter in Greek alphabet ) with footnote j = 1 and headnote k)(bj - aj ))

B = ( \Sigma( capitalized 18th letter in Greek alphabet) with footnote i = 1 and headhote n) f(a1 + (c((i - 1)*k + 1))*(b1 - a1),..., ak + (c(i*k))*(bk - ak)).
Yn = ( A/n )*B
the expense is O(n*k).
That�s it. And I think what we(my classmates and me) are really wildly guessing is the index of c in B. The three dots in B really causes headache. Still trying. Hopefully it can be done before the deadline.

Best Regards,
Ellen

[Jim Yingst]

I don't think the added dimensionality will change that a lot, except to make it smaller
No, each new dimension will increase the value of I by a factor of approximately pi/2, because for each new dimension you integrate along a new range of 0 to pi/2. So it's reasonable to expect for ten dimensions, the integral will be approximately (pi/2)^10 times as big as it was.
However, this ignores another effect - the added dimensions may also affect the average value of y in the region we're integrating over. Studying the function, we can see that each sine function will be somewhere between 0 and 1, and when we multiply ten sine functions we will get another value between 0 and 1 - but it's much more likely to be closer to 0 than 1, since if any of the 10 sines is 0 or near 0, the product of all ten will be at least that close to 0. So the exponential of the product of ten sines will be the exponential of a number which is probably very close to 0 - which means the exponential is probably very close to 1, slightly higher.
So - the expected value of this integral is
(range of x)^10 * (avg value of y)
(pi/2)^10 * 1
91.4524
This is pretty close to what the Monte Carlo integrals are converging to. It's a bit surprising to me that the average value of y should be that close to 1, since it does have a maximum value of exp(1) = e = 2.71. But I guess such high values are very rare in the range of this integral.

[Ellen Zhao]

Finally settled. Thanks to Mr. Yingst and Mr. Miller.

package Integration;
/**
 * Function integrator that implements universal high dimension
 * Monte Carlo algorithm.
 */
 
import java.util.Random;
public class HighDimensionMonteCarloIntegrator {
private static AlignRight ar = new AlignRight();
private static Random random = new Random();
    /**
     * Integrate the function from a to b using Monte Carlo algorithm.
     */
    public void integrate()
    {
    for ( int dotsNum = 100; dotsNum <= 10000000; dotsNum *= 10 )
    {
    double I         = 0;
    double sumf      = 0; // f is a auxiliary function
    double temp      = 1;
    int    dimension = 10;
        
    for (int j = 1; j <= dotsNum; ++j) {
    for (int i = 1; i <= dimension; ++i) {
    temp *= Math.sin( random.nextDouble()*Math.PI/2 );
    }    
    sumf += Math.exp( temp );
        
        if(j == dotsNum) {
        I = sumf*Math.pow(Math.PI/2, dimension)/j;
        ar.print(dotsNum, 9);
        ar.print(I, 25);
        ar.println();
        }         
    }
    }
 }
    
public static void main(String[] args) {
    ar.print("n", 9); ar.print("I", 25);
ar.underline();

HighDimensionMonteCarloIntegrator integrator =  
new HighDimensionMonteCarloIntegrator();
integrator.integrate();    
}
}  

The result is here:

        n                        I
----------------------------------
      100        91.45735204492146
     1000         91.4532756820613
    10000        91.45319580754229
   100000        91.45317456450654
  1000000        91.45317186857041
 10000000        91.45317140415479
Press any key to continue...  

Best Regards,
Ellen
[ January 18, 2003: Message edited by: Ellen Fu ]