• Post Reply Bookmark Topic Watch Topic
  • New Topic
programming forums Java Mobile Certification Databases Caching Books Engineering Micro Controllers OS Languages Paradigms IDEs Build Tools Frameworks Application Servers Open Source This Site Careers Other all forums
this forum made possible by our volunteer staff, including ...
Marshals:
  • Campbell Ritchie
  • Paul Clapham
  • Ron McLeod
  • Bear Bibeault
  • Liutauras Vilda
Sheriffs:
  • Jeanne Boyarsky
  • Tim Cooke
  • Junilu Lacar
Saloon Keepers:
  • Tim Moores
  • Tim Holloway
  • Stephan van Hulst
  • Jj Roberts
  • Carey Brown
Bartenders:
  • salvin francis
  • Frits Walraven
  • Piet Souris

High dimension Monte Carlo Integration

 
Ranch Hand
Posts: 581
  • Mark post as helpful
  • send pies
    Number of slices to send:
    Optional 'thank-you' note:
  • Quote
  • Report post to moderator
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 ]
 
Ranch Hand
Posts: 56
  • Mark post as helpful
  • send pies
    Number of slices to send:
    Optional 'thank-you' note:
  • Quote
  • Report post to moderator
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
Ranch Hand
Posts: 581
  • Mark post as helpful
  • send pies
    Number of slices to send:
    Optional 'thank-you' note:
  • Quote
  • Report post to moderator
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
Ranch Hand
Posts: 56
  • Mark post as helpful
  • send pies
    Number of slices to send:
    Optional 'thank-you' note:
  • Quote
  • Report post to moderator
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
Ranch Hand
Posts: 581
  • Mark post as helpful
  • send pies
    Number of slices to send:
    Optional 'thank-you' note:
  • Quote
  • Report post to moderator
Hi Norm,
Thank you so much! I�ll put the output here when I have completely implemented your idea.
Best Regards,
Ellen
 
Ellen Zhao
Ranch Hand
Posts: 581
  • Mark post as helpful
  • send pies
    Number of slices to send:
    Optional 'thank-you' note:
  • Quote
  • Report post to moderator
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
Ranch Hand
Posts: 56
  • Mark post as helpful
  • send pies
    Number of slices to send:
    Optional 'thank-you' note:
  • Quote
  • Report post to moderator
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
Ranch Hand
Posts: 581
  • Mark post as helpful
  • send pies
    Number of slices to send:
    Optional 'thank-you' note:
  • Quote
  • Report post to moderator
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
 
Wanderer
Posts: 18671
  • Mark post as helpful
  • send pies
    Number of slices to send:
    Optional 'thank-you' note:
  • Quote
  • Report post to moderator
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
Ranch Hand
Posts: 581
  • Mark post as helpful
  • send pies
    Number of slices to send:
    Optional 'thank-you' note:
  • Quote
  • Report post to moderator
Finally settled. Thanks to Mr. Yingst and Mr. Miller.

The result is here:

Best Regards,
Ellen
[ January 18, 2003: Message edited by: Ellen Fu ]
 
rubbery bacon. crispy tiny ad:
the value of filler advertising in 2020
https://coderanch.com/t/730886/filler-advertising
reply
    Bookmark Topic Watch Topic
  • New Topic