Monte Carlo Simulation
From Rosetta Code
Programming Task
This is a programming task. It lays out a problem which Rosetta Code users are encouraged to solve, using languages they know.
A simple Monte Carlo Simulation can be used to calculate the value for pi. If you had a circle and a square where the length of a side of the square was the same as the diameter of the circle, the ratio of the area of the circle to the area of the square would be pi/4. So, if you put this circle inside the square and select many random points inside the square, the number of points inside the circle divided by the number of points inside the square and the circle would be approximately pi/4.
Write a function to run a simulation like this with a variable number of random points to select. Also, show the results of a few different sample sizes.
Contents |
[edit] Ada
with Ada.Text_IO; use Ada.Text_IO; with Ada.Numerics.Float_Random; use Ada.Numerics.Float_Random; procedure Test_Monte_Carlo is Dice : Generator; function Pi (Throws : Positive) return Float is Inside : Natural := 0; begin for Throw in 1..Throws loop if Random (Dice) ** 2 + Random (Dice) ** 2 <= 1.0 then Inside := Inside + 1; end if; end loop; return 4.0 * Float (Inside) / Float (Throws); end Pi; begin Put_Line (" 10_000:" & Float'Image (Pi ( 10_000))); Put_Line (" 100_000:" & Float'Image (Pi ( 100_000))); Put_Line (" 1_000_000:" & Float'Image (Pi ( 1_000_000))); Put_Line (" 10_000_000:" & Float'Image (Pi ( 10_000_000))); Put_Line ("100_000_000:" & Float'Image (Pi (100_000_000))); end Test_Monte_Carlo;
The implementation uses built-in uniformly distributed on [0,1] random numbers. Note that the accuracy of the result depends on the quality of the pseudo random generator: its circle length and correlation to the function being simulated. Sample output:
10_000: 3.13920E+00
100_000: 3.14684E+00
1_000_000: 3.14197E+00
10_000_000: 3.14215E+00
100_000_000: 3.14151E+00
[edit] BASIC
Works with: QuickBasic version 4.5
Translation of: Java
DECLARE FUNCTION getPi! (throws!) CLS PRINT getPi(10000) PRINT getPi(100000) PRINT getPi(1000000) PRINT getPi(10000000) FUNCTION getPi (throws) inCircle = 0 FOR i = 1 TO throws 'a square with a side of length 2 centered at 0 has 'x and y range of -1 to 1 randX = (RND * 2) - 1'range -1 to 1 randY = (RND * 2) - 1'range -1 to 1 'distance from (0,0) = sqrt((x-0)^2+(y-0)^2) dist = SQR(randX ^ 2 + randY ^ 2) IF dist < 1 THEN 'circle with diameter of 2 has radius of 1 inCircle = inCircle + 1 END IF NEXT i getPi = 4! * inCircle / throws END FUNCTION
Output:
3.16 3.13648 3.142828 3.141679
[edit] D
import std.stdio: writefln; import std.random: rand; double pi(int nthrows) { int inside; for (int i; i < nthrows; i++) { double r1 = rand() / cast(double)uint.max; double r2 = rand() / cast(double)uint.max; if (r1*r1 + r2*r2 <= 1.0) inside++; } return 4.0 * inside / nthrows; } void main() { for (int n = 10_000; n <= 100_000_000; n *= 10) writefln("%9d: %07f", n, pi(n)); }
Sample output:
100000: 3.137760 1000000: 3.143508 10000000: 3.140836 100000000: 3.141666
For much faster results you can use the pseudorandom generator of Tango, or another one like R250-521.
[edit] Forth
Works with: GNU Forth
include random.fs 10000 value r : hit? ( -- ? ) r random dup * r random dup * + r dup * < ; : sims ( n -- hits ) 0 swap 0 do hit? if 1+ then loop ;
1000 sims 4 * . 3232 ok 10000 sims 4 * . 31448 ok 100000 sims 4 * . 313704 ok 1000000 sims 4 * . 3141224 ok 10000000 sims 4 * . 31409400 ok
[edit] Fortran
Works with: Fortran version 90 and later
MODULE Simulation
IMPLICIT NONE
CONTAINS
FUNCTION Pi(samples)
REAL :: Pi
REAL :: coords(2), length
INTEGER :: i, in_circle, samples
in_circle = 0
DO i=1, samples
CALL RANDOM_NUMBER(coords)
coords = coords * 2 - 1
length = SQRT(coords(1)*coords(1) + coords(2)*coords(2))
IF (length <= 1) in_circle = in_circle + 1
END DO
Pi = 4.0 * REAL(in_circle) / REAL(samples)
END FUNCTION Pi
END MODULE Simulation
PROGRAM MONTE_CARLO
USE Simulation
INTEGER :: n = 10000
DO WHILE (n <= 100000000)
WRITE (*,*) n, Pi(n)
n = n * 10
END DO
END PROGRAM MONTE_CARLO
Output:
10000 3.12120
100000 3.13772
1000000 3.13934
10000000 3.14114
100000000 3.14147
[edit] Haskell
import System.Random
import Control.Monad
get_pi throws = do results <- replicateM throws one_trial
return (4 * fromIntegral (sum results) / fromIntegral throws)
where
one_trial = do rand_x <- randomRIO (-1, 1)
rand_y <- randomRIO (-1, 1)
let dist = sqrt (rand_x*rand_x + rand_y*rand_y)
return (if dist < 1 then 1 else 0)
Example:
Prelude System.Random Control.Monad> get_pi 10000 3.1352 Prelude System.Random Control.Monad> get_pi 100000 3.15184 Prelude System.Random Control.Monad> get_pi 1000000 *** Exception: stack overflow
[edit] Java
public class MC { public static void main(String[] args) { System.out.println(getPi(10000)); System.out.println(getPi(100000)); System.out.println(getPi(1000000)); System.out.println(getPi(10000000)); System.out.println(getPi(100000000)); } public static double getPi(int numThrows){ int inCircle= 0; for(int i= 0;i < numThrows;i++){ //a square with a side of length 2 centered at 0 has //x and y range of -1 to 1 double randX= (Math.random() * 2) - 1;//range -1 to 1 double randY= (Math.random() * 2) - 1;//range -1 to 1 //distance from (0,0) = sqrt((x-0)^2+(y-0)^2) double dist= Math.sqrt(randX * randX + randY * randY); if(dist < 1){//circle with diameter of 2 has radius of 1 inCircle++; } } return 4.0 * inCircle / numThrows; } }
Output:
3.1396 3.14256 3.141516 3.1418692 3.14168604
[edit] MAXScript
fn monteCarlo iterations =
(
radius = 1.0
pointsInCircle = 0
for i in 1 to iterations do
(
testPoint = [(random -radius radius), (random -radius radius)]
if length testPoint <= radius then
(
pointsInCircle += 1
)
)
4.0 * pointsInCircle / iterations
)
[edit] OCaml
let get_pi throws = let rec helper i count = if i = throws then count else let rand_x = Random.float 2.0 -. 1.0 and rand_y = Random.float 2.0 -. 1.0 in let dist = sqrt (rand_x *. rand_x +. rand_y *. rand_y) in if dist < 1.0 then helper (i+1) (count+1) else helper (i+1) count in float (4 * helper 0 0) /. float throws
Example:
# get_pi 10000;; - : float = 3.15 # get_pi 100000;; - : float = 3.13272 # get_pi 1000000;; - : float = 3.143808 # get_pi 10000000;; - : float = 3.1421704 # get_pi 100000000;; - : float = 3.14153872
[edit] Python
[edit] At the interactive prompt
Python 2.6rc2 (r26rc2:66507, Sep 18 2008, 14:27:33) [MSC v.1500 32 bit (Intel)] on win32 IDLE 2.6rc2
>>> import random, math >>> throws = 1000 >>> 4.0 * sum(1 for p in xrange(throws) if math.hypot(*[random.random()*2-1 for q in [0,1]]) < 1) / float(throws) 3.1520000000000001 >>> throws = 1000000 >>> 4.0 * sum(1 for p in xrange(throws) if math.hypot(*[random.random()*2-1 for q in [0,1]]) < 1) / float(throws) 3.1396359999999999 >>> throws = 100000000 >>> 4.0 * sum(1 for p in xrange(throws) if math.hypot(*[random.random()*2-1 for q in [0,1]]) < 1) / float(throws) 3.1415666400000002
[edit] As a program using a function
from random import random from math import hypot try: import psyco psyco.full() except: pass def pi(nthrows): inside = 0 for i in xrange(nthrows): if hypot(random(), random()) < 1: inside += 1 return 4.0 * inside / nthrows for n in [10**4, 10**6, 10**7, 10**8]: print "%9d: %07f" % (n, pi(n))
Categories: Programming Tasks | Mathematical operations | Ada | BASIC | D | Forth | Fortran | Haskell | Java | MAXScript | OCaml | Python
