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.

Improve this question
8
  • 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. Commented Oct 11, 2015 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? Commented Oct 11, 2015 at 16:22
  • Remark: If you remove the square root, your program will run faster ;-) Commented Oct 11, 2015 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... Commented Oct 12, 2015 at 9:05
  • 2
    Every time someone teaches a novice to use real(8), god kills a kitten stackoverflow.com/questions/838310/fortran-90-kind-parameter Commented Oct 12, 2015 at 9:32

2 Answers 2

Reset to default
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).

Improve this answer
Sign up to request clarification or add additional context in comments.

1 Comment

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

Improve this answer

2 Comments

This works too, thank you. See the edit as my problem now is another one
have you tried precisions in between like 5E-3 or 3E-4, etc.?

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.