Monte Carlo pi method with fixed x

151 views Asked by At

In all Monte Carlo example codes I have found which calculate pi, both x and y are generated randomly between 0 and 1. For example, an example code looks like

  Ran rdm(time(NULL));
  double x, y;
  int sum=0;
  for(int i=0;i<N;i++)
  {
    x=rdm.doub();         // Both x and y are generated randomly
    y=rdm.doub();
    if(x*x+y*y<1)
      sum++;
  }
  return 4.*sum/N;

I don't understand what is the point to randomly generate both axis rather than choose a fixed and uniform x list and then only generate y randomly, like an example code below.

  Ran rdm(time(NULL));
  double x=0.5/N, dx=1./N, y;   //x is a fixed sequential list
  int sum=0;
  for(int i=0;i<N;i++)
  {
    y=rdm.doub();               // only y is generated randomly
    if(x*x+y*y<1)
      sum++;
    x+=dx;
  }
  return 4.*sum/N;

I tried both examples. Both methods converge with speed ~ 1/sqrt(N). The error of the second method is slightly smaller and the second code runs slightly faster. So it seems to me that the second method is better. And for sampling in higher dimension space, I think one can always pick a fixed list in one dimension and randomly sampling in other dimensions. Is this true?

Is there any reference that can prove the fixed x method is also correct?

If the fixed x method is better, why I have never seen any example code written in this way? Is it related to the adaptive importance sampling?

2

There are 2 answers

5
murrdpirate On BEST ANSWER

I believe it will work even if both x and y are generated uniformly. In that case, you'd be laying out a uniform grid of points and seeing which ones fall inside the circle.

I think you are right that for a given number of points, the uniform sampling would be as accurate as the random sampling, but significantly faster.

But what if you don't know the number of points you should use? In this case, I think the random sampling is better, because you can keep adding points until the overall value changes very little. You can't simply increment the number of points if you're using a uniform distribution. Any time you increase the number of points, you have to create an entire grid.

2
David Eisenstat On

Your method for computing π is indeed valid. I'm sure math.SE would give you a crisp proof of this, but it's basically (1) the unit disk is Riemann integrable (2) the indicator variable for each slice has the right expectation (3) Hoeffding's inequality (e.g.) to show convergence in probability for the sum.

Thing is, this method doesn't work on sets that are not Riemann integrable. Consider, for instance, the set ([0, 1] - Q) × ([0, 1] - Q), which has measure 1 but on which your method will always sample the rational x columns, which are not in the set.

See also low-discrepancy sequences, which is your idea taken a step further.