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 ]
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]
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
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
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.
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.
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.