Talk:Monte Carlo integration

merge VEGAS section with VEGAS algorithm

I suggest the entire section about VEGAS be merged with the VEGAS algorithm article. The two compliment each other and are similar in length. -Stimpy 16:23, 5 November 2006 (UTC)

Note also that the VEGAS section of this article is copied verbatim from the GSL help pages about the algorithm: [1]. DutchCanadian (talk) 21:58, 13 March 2013 (UTC)

Also it is disproportionately long in this article compared to the importance of the VEGAS algorithm. Livingthingdan (talk) 07:23, 8 April 2020 (UTC)

Parameters, etc.

This article should focus on the general principles of each method, not the documentation of a specific implementation. Names of parameters in a software package have no place in a general article about Monte Carlo integration, nor anywhere else in Wikipedia. This is an encyclopedia, not a documentation archive. Mtford 23:28, 5 October 2007 (UTC)

You're of course right. The problem is that nobody has yet rewritten the article. It would be great if you could do this. -- Jitse Niesen (talk) 01:10, 6 October 2007 (UTC)

made a commit on this line. I would like to confirm that nothing really important was erased. I tried to remove the details of the implementation, and keep the big ideas and lines of though of the method in here. Jorgecarleitao (talk) 16:49, 22 July 2012 (UTC)

pi/4 = .785?

• The area of the square is 4
• The area of the circle is pi*(1*1)
• The area of the circle is 80% of the area of the square, they say.
• AreaCircle = pi = 80%*4 = 3.2
• pi/4 = 3.2/4 = .8

How do they get pi/4 = .785? Martin | talk • contribs 05:07, 22 July 2009 (UTC)

pi = (pi/4) * 4 = 0.785 * 4 = 3.14 lol —Preceding unsigned comment added by 95.181.12.52 (talk) 21:24, 16 December 2009 (UTC)

• Furthermore, what is the point of using monte carlo to approximate pi if you have to use pi in the calculation? —Preceding unsigned comment added by 20.132.68.146 (talk) 19:04, 23 March 2010 (UTC)

The approximation that would be achieved to pi/4 would be 0.8, as is mentioned above. It would take more sample points to achieve the approximation that is in that caption, I will fix that now. The calculation of the value of pi using monte carlo integration does not require using pi in the calculation at all. jgoldfar (talk) 15:19, 22 December 2010 (UTC)

JavaScript example

A chunk of code has been removed from the article. I think it could be useful if it was re-written to be more readable, as this integration algorithm is used strictly in numerical applications due to number of steps required. http://en.wiki.x.io/w/index.php?title=Monte_Carlo_integration&action=historysubmit&diff=337812469&oldid=337197640

131.111.8.104 (talk) 05:26, 15 January 2010 (UTC)

Difference between Monte Carlo integration and Monte Carlo methods unclear

"Monte Carlo integration is numerical integration using random numbers." ... "Monte Carlo methods, however, randomly choose the points at which the integrand is evaluated." Isn't that the same thing? Or am I missing something? —Kri (talk) 21:48, 14 July 2011 (UTC)

changed introductory text to be more specific on what it is supposed to say. Jorgecarleitao (talk) 16:51, 22 July 2012 (UTC)

Two different algorithms?

I think this article actually states two different algorithms for plain (I think "naive" is the word that's usually used?) monte carlo integration. The first is what's specified in the introduction: "Informally, to estimate the area of a domain D, first pick a simple domain E whose area is easily calculated and which contains D. Now pick a sequence of random points that fall within E. Some fraction of these points will also fall within D. The area of D is then estimated as this fraction multiplied by the area of E." This first algorithm corresponds to the picture/caption of the circle.

In other words, that algorithm is this: Define a hypercube that encloses your entire integral. Pick points randomly within that hypercube. Determine the average number of those points which fall within the surface defined by your multi-dimensional function, and multiply that by the volume (hypervolume?) of the hypercube to get the integral. This involves randomly picking points in the output dimensions, and comparing them to results of the function when evaluated on just the input dimensions.

However, the algorithm that's actually specified in section 1 is this: Define a hypercube that enclose only the region of the integral (not the inputs). Pick values randomly for each input dimension (but not output dimension). Evaluate your multidimensional function at each of these points, average them, then multiply them by the volume of the hypercube region.

Let's say you have a function y = f(x), and you want to compute the integral of f(x) with respect to x, over the range a to b.

The first algorithm would require defining a rectangle with the left edge of the rectangle at x = a, the right edge at x = b, the bottom at y = i, and top at y = h, where you know that i <= f(x) <= h for all values of x between a and b. Then you randomly pick points (x,y) within the rectangle. You find the probability of a random point being within the function (i.e. for a random point (x,y) within the rectangle you find P(f(x) >= y)), and multiply that by the area of the rectangle.

The second algorithm would simply randomly pick values for x between a and b, find the average value for f(x) of the points, and multiply that by b - a.

Although both algorithms will converge on a correct answer, the first one requires being able to set an upper bound on a function within a range, and it is generally far less efficient. If the first algorithm isn't removed from the article altogether, the distinction between them should be made clear.

Novakog (talk) 08:35, 24 July 2011 (UTC)

Mathematica example

I agree with the addition of mathematica code to this page, to show an example of implementation. However, I think no code describes a method: the method is described by a mathematical formula/procedure, which has to be independent of the programming language to be used.

This said, the text before the actual code must explain 1. what one wants to exemplify; 2. what the code does. Every line of code should also have comments (*comment*), every time is not explicit what it means.Jorgecarleitao (talk) 08:19, 19 December 2012 (UTC)

Made a cleanup on the code, which had some minor errors, and suggest to change it for the following:

(* function to be integrated, on interval [0.8,3] *) func[x_] := 1/(1 + Sinh[2 x] Log[x]^2); p = Plot[func[x], {x, 0.8, 3}, PlotStyle -> Black]; p1 = Plot[PDF[NormalDistribution[1, 0.4], x], {x, 0.8, 3}, PlotStyle -> Blue, Filling -> Bottom]; p2 = Plot[PDF[UniformDistribution[{0.8, 3}], x], {x, 0.8, 3}, PlotStyle -> Red, Filling -> Bottom]; Show[{p, p1, p2}] NSolve[D[func[x], x] == 0, x, Reals]; (* the maximum of the function is at 1 *) (*** parameters for the random generation ***) (* use n points \ generated from a normal distribution with \[Mu]=1, \[Sigma]=0.4 \ truncated on interval [0.8,3] *) n = 1000; normal["points"] = RandomVariate[ TruncatedDistribution[{0.8, 3}, NormalDistribution[1, 0.4]], n]; normal["normalization"] = Integrate[PDF[NormalDistribution[1, 0.4], x], {x, 0.8, 3}]; (* use n points generated from a uniform distribution on interval \ [0.8,3] *) uniform["points"] = RandomVariate[UniformDistribution[{0.8, 3}], n]; uniform["normalization"] = 3 - 0.8; (* results: expected, by normal distribution, by uniform distribution*) \ NIntegrate[func[x], {x, 0.8, 3}] normal["integral"] = Total[func[ normal["points"]]/(PDF[NormalDistribution[1, 0.4], normal["points"]]/normal["normalization"])]/n uniform["integral"] = Total[func[uniform["points"]]/(1/uniform["normalization"])]/n

However, I think mathematica code on wikipedia, as it is, is really ugly. There should be a way of having it with a cleaner presentation.Jorgecarleitao (talk) 09:14, 19 December 2012 (UTC)

I just ran the example in R and had to scale by 8 instead of 4 to get a more accurate result

I am not sure why, but the example with a 4 just didn't work out; is there a scaling issue here? Here is my R code

thefunc <- function(point){
  output = point[1]^2+point[2]^2
  if(output<=1){
    output
  }
  else 0
}

N = 100000
samplepoints = matrix(nrow = N,ncol=2,data=runif(2*N,-1,1))
estimate = 8*mean(apply(samplepoints,1,thefunc))

EDIT: I'm an idiot and figured this out. The function in my code was wrong.

thefunc <- function(point){
  if(point[1]^2+point[2]^2<=1) 1
  else 0
}

Overview section describes a very specific case

The overview section describes a very specific case, where (a) the integration space is finite (b) the sampling is unbiased ("naive"). This hides an important general concept that belongs in the overview, that the weight of the i-th event (drawn at x_i) is 1/PDF(x_i).

It could be fine to start a lesson with such a simple case. But it's should be introduced as such. On the other hand, starting an encyclopedic entry with an overview is a fine idea, and giving the general case in a similar style seems like it would work well.