4

I'm taking a course in Numerical Methods and I've been requested to implement the famous Monte Carlo algorithm to find pi that you can find here.

I had no difficulties in writing the code with an arbitrary number of trials:

REAL(8) FUNCTION distance(xvalue, yvalue) RESULT(dist)
IMPLICIT NONE
REAL(8), INTENT(in) :: xvalue, yvalue
dist = SQRT(xvalue**2 + yvalue**2)
END FUNCTION distance 

PROGRAM ass2
  IMPLICIT NONE

  INTEGER, DIMENSION(1) :: SEED
  REAL(8) :: p, x, y

  REAL(8), EXTERNAL :: distance

  REAL(8) :: pi_last, pi
  INTEGER :: npc, npt, i

  npc = 0
  npt = 0
  pi = 1.0

  SEED(1) = 12345
  CALL RANDOM_SEED

  DO i=1, 1000000000

     CALL RANDOM_NUMBER(p)
     x = p
     CALL RANDOM_NUMBER(p)
     y = p

     npt = npt + 1

     IF (distance(x, y) < 1.0)  THEN
          npc = npc + 1
     END IF

     pi_last = pi
     pi = 4.0*(npc*1.0)/(npt*1.0)

  END DO

  PRINT*, 'Pi:', pi

END PROGRAM ass2

I noticed that it converges approximately as sqrt(N of steps). Now I have to stop the algorithm at a certain precision, so I created an endless DO loop with an EXIT inside an IF statement:

    REAL(8) FUNCTION distance(xvalue, yvalue) RESULT(dist)
IMPLICIT NONE
REAL(8), INTENT(in) :: xvalue, yvalue
dist = SQRT(xvalue**2 + yvalue**2)
END FUNCTION distance 

PROGRAM ass2
  IMPLICIT NONE

  INTEGER, DIMENSION(1) :: SEED
  REAL(8) :: p, x, y

  REAL(8), EXTERNAL :: distance

  REAL(8) :: pi_last, pi
  INTEGER :: npc, npt, i

  npc = 0
  npt = 0
  pi = 1.0

  SEED(1) = 12345
  CALL RANDOM_SEED

  DO

     CALL RANDOM_NUMBER(p)
     x = p
     CALL RANDOM_NUMBER(p)
     y = p

     npt = npt + 1

     IF (distance(x, y) < 1.0)  THEN
          npc = npc + 1
     END IF

     pi_last = pi
     pi = 4.0*(npc*1.0)/(npt*1.0)

  IF ( ABS(pi - pi_last) < 0.000001 .AND. pi - pi_last /= 0)    THEN
    EXIT
  END IF

  END DO

  PRINT*, 'Pi:', pi

END PROGRAM ass2

The problem is that this returns a value of pi that doesn't have the precision I asked for. I get the logic behind it: if I get two consecutive values far from pi but close to each other, the condition will be satisfied and the program will exit the DO statement. The problem is that I don't know how to modify it in order to get a precision decided by me. So the question is:

How do I implement this algorithm in a manner such that I can decide the precision of pi in output?

EDIT: Ok, I implemented both of your solutions and they work, but only for 10^(-1), 10^(-3) and 10^(-5). I think it's a problem of the pseudorandom sequence if for 10^(-2) and 10^(-4) it returns an incorrect value of pi.

iacolippo
  • 4,133
  • 25
  • 37
  • Without going through your code, your precision is hardcoded to be 0.000001, does it achieve that? Note, that this measure is not in relation to Pi, but two successive iterations. If you want this to be changeable at runtime, you need a variable and get the value either as an argument to your command or read it in the initialization. Look for "read" for the last and "get_command_argument" for the first variant. – haraldkl Oct 11 '15 at 14:08
  • I think I've not been clear, sorry, my fault. The problem is that it doesn't achieve that hardcoded precision. I don't want to be able to change the value at runtime, neither it's necessary that I get it from keyboard with READ* command :) I noticed too that it's in relation to two successive iterations, but note that I can't define a constant pi with value from books or web to check if the digit is correct. I have to get pi with a certain precision without knowing its value. How do I do it? – iacolippo Oct 11 '15 at 16:22
  • Remark: If you remove the square root, your program will run faster ;-) – Anthony Scemama Oct 11 '15 at 19:23
  • One usually monitors the convergence of MC by calculating statistical error somehow. Because the samples here are statistically independent, the standard deviation can be used for this purpose. Also, it should be better to use double-precision literals like 1.0d0... – roygvib Oct 12 '15 at 09:05
  • 2
    Every time someone teaches a novice to use `real(8)`, god kills a kitten http://stackoverflow.com/questions/838310/fortran-90-kind-parameter – Vladimir F Героям слава Oct 12 '15 at 09:32
  • LOL I will use selected kind starting by now – iacolippo Oct 12 '15 at 09:38
  • all else aside, your convergence criteria isn't a very good one (looking only at the last two points). In this case since you know the answer why don't you check against a high precision value for `Pi` ? If you did not know the result you would need a more sophisticated criteria, tracking the last few hundred results for example. – agentp Oct 13 '15 at 16:01
  • Yes, I actually wrote in the question that my criteria wasn't good :) Even the solution you suggest (of tracking the last few hundred) doesn't work, and I should have some criteria on how to choose this number of steps to check. The other solution, of checking against a high precision value, isn't allowed, I have to do this without knowing the value of pi. I've delivered my homework, I will post here if I get the solution. Thank to you all :) – iacolippo Oct 14 '15 at 16:40

2 Answers2

4

It's not enough to specify a desired precision -- you also need to allow some chance that the precision target is not met. Then you can solve for the number of trials in (e.g.) Hoeffding's inequality to meet your desired precision with the desired probability (as you've observed, n needs to be about the square root of one over the precision to succeed with constant probability).

David Eisenstat
  • 64,237
  • 7
  • 60
  • 120
  • This makes sense :) It works for 10^(-3) and 10^(-5), but not for 10^(-4) and I think it is a problem of the pseudorandom number generator that gives a sequence that falls in the 0.1% of probability the target isn't met. I'm new to FORTRAN so I just took it as the teacher gave it to me. – iacolippo Oct 11 '15 at 18:57
3

In an ideal setting (perfect random number generator generating real numbers in the mathematical sense) your variable npc is a random variable with binomial distribution B(n,π/4) where n is npt from your program. Its expected value is n*π/4, so you correctly compute the approximation of π as pi=4*npc/npt. Now this approximation can take all values from 0 to 4 no matter how many loop iterations you calculate, because npc can take all values from 0 to npt. For a range around π you can only give a probability (using c as shorthand for npc; P denotes the probability of an event):

P(|pi - π| < d) = P(-d < pi - π < d) = P(-d < 4*c/n - π < d) = P(n*(π-d)/4 < c < n*(π+d)/4) =

= P(c < n*(π+d)/4) - P(c < n*(π-d)/4) ~=

~= FN(n*(π+d)/4) - FN(n*(π-d)/4) = 2F(d*√(n/(π(4-π))))-1

where FN is the probability function of the normal distribution N(nπ/4;nπ/4(1-π/4)) which approximates the above binomial distribution and F is the probability function of the standard normal distribution. Now given a deviation d and a probability p you can compute n s.t. the last term does not go below p:

n = ceil(π(4-π)(F-1((p+1)/2)/d)^2))

Then with n loop iterations you can compute the approximation pi of π to a desired precision with given probability. If we want to achieve a probability of p=99%, then the above formula simplifies to

n ~= 17.89/d2 ,

so for precision d=0.0001 roughly n=1.789E9 iterations are necessary!

Note: since a computer cannot meet the above ideal setting there is also a (theoretical) limit on the precision you can reach with this algorithm. There are only finitely many floating point numbers representable in a computer, so your points (x,y) lie on a kind of grid. The best approximation of π that can be computed with this algorithm boils down to performing your loop over all grid points in [0,1]x[0,1]. The good old C-function rand() has a resolution of 31 bits (at least in the VS stdlib). So it does not make sense to compute more than n=312 points, which gives a maximal precision of √(17.89/n) = 1.97E-9 when asking for 99% correctness.

coproc
  • 6,027
  • 2
  • 20
  • 31