Monte Carlo methods
From Rosetta Code
You are encouraged to solve this task according to the task description, using any language you may know.
A simple Monte Carlo Simulation can be used to calculate the value for π. 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 π/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 π/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. For software where the number π is not built-in, we give π to a couple of digits: 3.141592653589793238462643383280
[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] ALGOL 68
PROC pi = (INT throws)REAL:
BEGIN
INT inside := 0;
TO throws DO
IF random ** 2 + random ** 2 <= 1 THEN
inside +:= 1
FI
OD;
4 * inside / throws
END # pi #;
print ((" 10 000:",pi ( 10 000),new line));
print ((" 100 000:",pi ( 100 000),new line));
print ((" 1 000 000:",pi ( 1 000 000),new line));
print ((" 10 000 000:",pi ( 10 000 000),new line));
print (("100 000 000:",pi (100 000 000),new line))
Sample output:
10 000:+3.15480000000000e +0
100 000:+3.12948000000000e +0
1 000 000:+3.14169200000000e +0
10 000 000:+3.14142040000000e +0
100 000 000:+3.14153276000000e +0
[edit] AutoHotkey
Search autohotkey.com: Carlo methods
Source: AutoHotkey forum by Laszlo
MsgBox % MontePi(10000) ; 3.154400
MsgBox % MontePi(100000) ; 3.142040
MsgBox % MontePi(1000000) ; 3.142096
MontePi(n) {
Loop %n% {
Random x, -1, 1.0
Random y, -1, 1.0
p += x*x+y*y < 1
}
Return 4*p/n
}
[edit] AWK
# --- with command line argument "throws" ---
BEGIN{ th=ARGV[1];
for(i=0; i<th; i++) cin += (rand()**2 + rand()**2) < 1 ;
printf("Pi = %8.5f\n",4*cin/th);
}
usage: awk -f pi 2300
Pi = 3.14333
[edit] BASIC
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] BBC BASIC
PRINT FNmontecarlo(1000)
PRINT FNmontecarlo(10000)
PRINT FNmontecarlo(100000)
PRINT FNmontecarlo(1000000)
PRINT FNmontecarlo(10000000)
END
DEF FNmontecarlo(t%)
LOCAL i%, n%
FOR i% = 1 TO t%
IF RND(1)^2 + RND(1)^2 < 1 n% += 1
NEXT
= 4 * n% / t%
Sample output:
3.136
3.1396
3.13756
3.143624
3.1412816
[edit] C
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
double pi(double tolerance)
{
double x, y, val, error;
unsigned long sampled = 0, hit = 0, i;
do {
/* don't check error every turn, make loop tight */
for (i = 1000000; i; i--, sampled++) {
x = rand() / (RAND_MAX + 1.0);
y = rand() / (RAND_MAX + 1.0);
if (x * x + y * y < 1) hit ++;
}
val = (double) hit / sampled;
error = val * sqrt(val * (1 - val) / sampled) * 4;
val *= 4;
/* some feedback, or user gets bored */
fprintf(stderr, "Pi = %f +/- %5.3e at %ldM samples.\r",
val, error, sampled/1000000);
} while (!hit || error > tolerance);
/* !hit is for completeness's sake; if no hit after 1M samples,
your rand() is BROKEN */
return val;
}
int main()
{
printf("Pi is %f\n", pi(3e-4)); /* set to 1e-4 for some fun */
return 0;
}
[edit] C#
using System;
class Program {
static double MonteCarloPi(int n) {
int inside = 0;
Random r = new Random();
for (int i = 0; i < n; i++) {
if (Math.Pow(r.NextDouble(), 2)+ Math.Pow(r.NextDouble(), 2) <= 1) {
inside++;
}
}
return 4.0 * inside / n;
}
static void Main(string[] args) {
int value = 1000;
for (int n = 0; n < 5; n++) {
value *= 10;
Console.WriteLine("{0}:{1}", value.ToString("#,###").PadLeft(11, ' '), MonteCarloPi(value));
}
}
}
Output
10,000:3.1436
100,000:3.14632
1,000,000:3.139476
10,000,000:3.1424476
100,000,000:3.1413976
[edit] Clojure
(defn calc-pi [iterations]
(loop [x (rand) y (rand) in 0 total 1]
(if (< total iterations)
(recur (rand) (rand) (if (<= (+ (* x x) (* y y)) 1) (inc in) in) (inc total))
(double (* (/ in total) 4)))))
(doseq [x (take 5 (iterate #(* 10 %) 10))] (println (str (format "% 8d" x) ": " (calc-pi x))))
output:
100: 3.2 1000: 3.124 10000: 3.1376 100000: 3.14104 1000000: 3.141064
[edit] Common Lisp
(defun approximate-pi (n)
(/ (loop repeat n count (<= (abs (complex (random 1.0) (random 1.0))) 1.0)) n 0.25))
(dolist (n (loop repeat 5 for n = 1000 then (* n 10) collect n))
(format t "~%~8d -> ~f" n (approximate-pi n)))
Output:
1000 -> 3.132
10000 -> 3.1184
100000 -> 3.1352
1000000 -> 3.142072
10000000 -> 3.1420677
[edit] D
import std.stdio, std.random, std.math;
double pi(in int nthrows) {
int inside;
foreach (i; 0 .. nthrows)
if (hypot(uniform(0, 1.0), uniform(0, 1.0)) <= 1)
inside++;
return 4.0 * inside / nthrows;
}
void main() {
foreach (p; 1 .. 8)
writefln("%10s: %07f", 10 ^^ p, pi(10 ^^ p));
}
Sample output:
10: 3.200000
100: 3.120000
1000: 3.076000
10000: 3.140400
100000: 3.146520
1000000: 3.140192
10000000: 3.141476
Sample output with foreach(p;1..10):
10: 3.200000
100: 3.240000
1000: 3.180000
10000: 3.150400
100000: 3.143080
1000000: 3.140996
10000000: 3.141442
100000000: 3.141439
1000000000: 3.141559
[edit] More Functional Style
Similar output.
import std.stdio, std.random, std.math, std.algorithm, std.range;
bool isIn(int){ return hypot(uniform(0,1.0), uniform(0,1.0)) <= 1; }
double pi(in int n) { return 4.0 * count!isIn(iota(n)) / n; }
void main() {
foreach (p; 1 .. 8)
writefln("%10s: %07f", 10 ^^ p, pi(10 ^^ p));
}
[edit] E
This computes a single quadrant of the described square and circle; the effect should be the same since the other three are symmetric.
def pi(n) {
var inside := 0
for _ ? (entropy.nextFloat() ** 2 + entropy.nextFloat() ** 2 < 1) in 1..n {
inside += 1
}
return inside * 4 / n
}
Some sample runs:
? pi(10) # value: 2.8 ? pi(10) # value: 2.0 ? pi(100) # value: 2.96 ? pi(10000) # value: 3.1216 ? pi(100000) # value: 3.13088
? pi(100000) # value: 3.13848
[edit] Euler Math Toolbox
>function map MonteCarloPI (n,plot=false) ...
$ X:=random(1,n);
$ Y:=random(1,n);
$ if plot then
$ plot2d(X,Y,>points,style=".");
$ plot2d("sqrt(1-x^2)",color=2,>add);
$ endif
$ return sum(X^2+Y^2<1)/n*4;
$endfunction
>MonteCarloPI(10^(1:7))
[ 3.6 2.96 3.224 3.1404 3.1398 3.141548 3.1421492 ]
>pi
3.14159265359
>MonteCarloPI(10000,true):
[edit] Factor
Since Factor lets the user choose the range of the random generator, we use 2^32.
USING: kernel math math.functions random sequences ;
: limit ( -- n ) 2 32 ^ ; inline
: in-circle ( x y -- ? ) limit [ sq ] tri@ [ + ] [ <= ] bi* ;
: rand ( -- r ) limit random ;
: pi ( n -- pi ) [ [ drop rand rand in-circle ] count ] keep / 4 * >float ;
Example use:
10000 pi .
3.1412
[edit] Fantom
class MontyCarlo
{
// assume square/circle of width 1 unit
static Float findPi (Int samples)
{
Int insideCircle := 0
samples.times
{
x := Float.random
y := Float.random
if ((x*x + y*y).sqrt <= 1.0f) insideCircle += 1
}
return insideCircle * 4.0f / samples
}
public static Void main ()
{
[100, 1000, 10000, 1000000, 10000000].each |sample|
{
echo ("Sample size $sample gives PI as ${findPi(sample)}")
}
}
}
Output:
Sample size 100 gives PI as 3.2 Sample size 1000 gives PI as 3.132 Sample size 10000 gives PI as 3.1612 Sample size 1000000 gives PI as 3.139316 Sample size 10000000 gives PI as 3.1409272
[edit] 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
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] Go
package main
import (
"fmt"
"math"
"math/rand"
"time"
)
func getPi(numThrows int) float64 {
inCircle := 0
for 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
randX := rand.Float64()*2 - 1 //range -1 to 1
randY := rand.Float64()*2 - 1 //range -1 to 1
//distance from (0,0) = sqrt((x-0)^2+(y-0)^2)
dist := math.Hypot(randX, randY)
if dist < 1 { //circle with diameter of 2 has radius of 1
inCircle++
}
}
return 4 * float64(inCircle) / float64(numThrows)
}
func main() {
rand.Seed(time.Now().UnixNano())
fmt.Println(getPi(10000))
fmt.Println(getPi(100000))
fmt.Println(getPi(1000000))
fmt.Println(getPi(10000000))
fmt.Println(getPi(100000000))
}
Output:
3.1164 3.1462 3.142892 3.1419692 3.14149596
[edit] Haskell
import System.Random
import Control.Monad
get_pi throws = do results <- replicateM throws one_trial
return (4 * fromIntegral (foldl' (+) 0 results) / fromIntegral throws)
where
one_trial = do rand_x <- randomRIO (-1, 1)
rand_y <- randomRIO (-1, 1)
let dist :: Double
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 3.145024
[edit] HicEst
FUNCTION Pi(samples)
inside = 0
DO i = 1, samples
inside = inside + ( (RAN(1)^2 + RAN(1)^2)^0.5 <= 1)
ENDDO
Pi = 4 * inside / samples
END
WRITE(ClipBoard) Pi(1E4) ! 3.1504
WRITE(ClipBoard) Pi(1E5) ! 3.14204
WRITE(ClipBoard) Pi(1E6) ! 3.141672
WRITE(ClipBoard) Pi(1E7) ! 3.1412856
[edit] Icon and Unicon
procedure main()Output:
every t := 10 ^ ( 5 to 9 ) do
printf("Rounds=%d Pi ~ %r\n",t,getPi(t))
end
link printf
procedure getPi(rounds)
incircle := 0.
every 1 to rounds do
if 1 > sqrt((?0 * 2 - 1) ^ 2 + (?0 * 2 - 1) ^ 2) then
incircle +:= 1
return 4 * incircle / rounds
end
Rounds=100000 Pi ~ 3.143400 Rounds=1000000 Pi ~ 3.141656 Rounds=10000000 Pi ~ 3.140437 Rounds=100000000 Pi ~ 3.141375 Rounds=1000000000 Pi ~ 3.141604
[edit] J
Explicit Solution:
piMC=: monad define "0
4* y%~ +/ 1>: %:+/*: <:+: (2,y) ?@$ 0
)
Tacit Solution:
piMCt=: (0.25&* %~ +/@(1 >: [: +/&.:*: _1 2 p. 0 ?@$~ 2&,))"0
Examples:
piMC 1e6
3.1426
piMC 10^i.7
4 2.8 3.24 3.168 3.1432 3.14256 3.14014
[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);
//^ or in Java 1.5+: double dist= Math.hypot(randX, 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] JavaScript
function mcpi(n){
var x,y,m=0;
for(var i = 0; i < n; i += 1) {
x = Math.random();
y = Math.random();
if (x*x + y*y < 1) { m += 1; }
}
return 4*m/n;
}
console.log(mcpi(1000));
console.log(mcpi(10000));
console.log(mcpi(100000));
console.log(mcpi(1000000));
console.log(mcpi(10000000));
3.168 3.1396 3.13692 3.140512 3.1417656
[edit] K
sim:{4*(+/{~1<+/(2_draw 0)^2}'!x)%x}
sim 10000
3.103
sim'10^!8
4 2.8 3.4 3.072 3.1212 3.14104 3.14366 3.1413
[edit] Liberty BASIC
for pow = 2 to 6
n = 10^pow
print n, getPi(n)
next
end
function getPi(n)
incircle = 0
for throws=0 to n
scan
incircle = incircle + (rnd(1)^2+rnd(1)^2 < 1)
next
getPi = 4*incircle/throws
end function
Sample output: 100 2.89108911 1000 3.12887113 10000 3.13928607 100000 3.13864861 1000000 3.13945686
[edit] Locomotive Basic
10 mode 1:randomize time:defint a-z
20 input "How many samples";n
30 u=n/100+1
40 r=100
50 for i=1 to n
60 if i mod u=0 then locate 1,3:print using "##% done"; i/n*100
70 x=rnd*2*r-r
80 y=rnd*2*r-r
90 if sqr(x*x+y*y)<r then m=m+1
100 next
110 pi2!=4*m/n
120 locate 1,3
130 print m;"points in circle"
140 print "Computed value of pi:"pi2!
150 print "Difference to real value of pi: ";
160 print using "+#.##%"; (pi2!-pi)/pi*100
[edit] Logo
to square :n
output :n * :n
end
to trial :r
output less? sum square random :r square random :r square :r
end
to sim :n :r
make "hits 0
repeat :n [if trial :r [make "hits :hits + 1]]
output 4 * :hits / :n
end
show sim 1000 10000 ; 3.18
show sim 10000 10000 ; 3.1612
show sim 100000 10000 ; 3.145
show sim 1000000 10000 ; 3.140828
[edit] LSL
To test it yourself; rez a box on the ground, and add the following as a New Script. (Be prepared to wait... LSL can be slow, but the Servers are typically running thousands of scripts in parallel so what do you expect?)
integer iMIN_SAMPLE_POWER = 0;
integer iMAX_SAMPLE_POWER = 6;
default {
state_entry() {
llOwnerSay("Estimating Pi ("+(string)PI+")");
integer iSample = 0;
for(iSample=iMIN_SAMPLE_POWER ; iSample<=iMAX_SAMPLE_POWER ; iSample++) {
integer iInCircle = 0;
integer x = 0;
integer iMaxSamples = (integer)llPow(10, iSample);
for(x=0 ; x<iMaxSamples ; x++) {
if(llSqrt(llPow(llFrand(2.0)-1.0, 2.0)+llPow(llFrand(2.0)-1.0, 2.0))<1.0) {
iInCircle++;
}
}
float fPi = ((4.0*iInCircle)/llPow(10, iSample));
float fError = llFabs(100.0*(PI-fPi)/PI);
llOwnerSay((string)iSample+": "+(string)iMaxSamples+" = "+(string)fPi+", Error = "+(string)fError+"%");
}
llOwnerSay("Done.");
}
}
Output:
Estimating Pi (3.141593) 0: 1 = 4.000000, Error = 27.323954% 1: 10 = 4.000000, Error = 27.323954% 2: 100 = 2.880000, Error = 8.326753% 3: 1000 = 3.188000, Error = 1.477192% 4: 10000 = 3.133600, Error = 0.254414% 5: 100000 = 3.138840, Error = 0.087620% 6: 1000000 = 3.142684, Error = 0.034739% Done.
[edit] Lua
function MonteCarlo ( n_throws )
math.randomseed( os.time() )
n_inside = 0
for i = 1, n_throws do
if math.random()^2 + math.random()^2 <= 1.0 then
n_inside = n_inside + 1
end
end
return 4 * n_inside / n_throws
end
print( MonteCarlo( 10000 ) )
print( MonteCarlo( 100000 ) )
print( MonteCarlo( 1000000 ) )
print( MonteCarlo( 10000000 ) )
Output:
3.1436 3.13636 3.14376 3.1420188
[edit] Mathematica
We define a function with variable sample size:
MonteCarloPi[samplesize_Integer] := N[4Mean[If[# > 1, 0, 1] & /@ Norm /@ RandomReal[1, {samplesize, 2}]]]
Example (samplesize=10,100,1000,....10000000):
{#, MonteCarloPi[#]} & /@ (10^Range[1, 7]) // Grid
gives back:
10 3.2
100 3.16
1000 3.152
10000 3.1228
100000 3.14872
1000000 3.1408
10000000 3.14134
monteCarloPi = 4. Mean[UnitStep[1 - Total[RandomReal[1, {2, #}]^2]]] &;
monteCarloPi /@ (10^Range@6)
[edit] MATLAB
See: Monte Carlo Simulation in MATLAB for more examples
The first example given is not vectorized. MATLAB has a self-imposed memory limitation that prevents this simulation from having more than 3 decimal digits of accuracy. Because of this limitation it is best to vectorize the code as much as possible so extra memory isn't consumed by unneeded variables. Therefore, I have provided a second solution that is maximally vectorized.
Minimally Vectorized:
function piEstimate = monteCarloPi(numDarts)
%The square has a sides of length 2, which means the circle has radius
%1.
%Generate a table of random x-y value pairs in the range [0,1] sampled
%from the uniform distribution for each axis.
darts = rand(numDarts,2);
%Any darts that are in the circle will have position vector whose
%length is less than or equal to 1 squared.
dartsInside = ( sum(darts.^2,2) <= 1 );
piEstimate = 4*sum(dartsInside)/numDarts;
end
Completely Vectorized:
function piEstimate = monteCarloPi(numDarts)
piEstimate = 4*sum( sum(rand(numDarts,2).^2,2) <= 1 )/numDarts;
end
Sample Output:
>> monteCarloPi(7000000)
ans =
3.141512000000000
[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] Octave
function p = montepi(samples)
in_circle = 0;
for samp = 1:samples
v = [ unifrnd(-1,1), unifrnd(-1,1) ];
if ( v*v.' <= 1.0 )
in_circle++;
endif
endfor
p = 4*in_circle/samples;
endfunction
l = 1e4;
while (l < 1e7)
disp(montepi(l));
l *= 10;
endwhile
Since it runs slow, I've stopped it at the second iteration, obtaining:
3.1560 3.1496
[edit] Much faster implementation
function result = montepi(n)
result = sum(rand(1,n).^2+rand(1,n).^2<1)/n*4;
endfunction
[edit] PARI/GP
MonteCarloPi(tests)=4.*sum(i=1,tests,norml2([random(1.),random(1.)])<1)/tests;
A hundred million tests (about a minute) yielded 3.14149000, slightly more precise (and round!) than would have been expected. A million gave 3.14162000 and a thousand 3.14800000.
[edit] Pascal
Program MonteCarlo(output);
uses
Math;
function MC_Pi(expo: integer): real;
var
x, y: real;
i, hits, samples: longint;
begin
samples := 10**expo;
hits := 0;
randomize;
for i := 1 to samples do
begin
x := random;
y := random;
if sqrt(x*x + y*y) < 1.0 then
inc(hits);
end;
MC_Pi := 4.0 * hits / samples;
end;
var
i: integer;
begin
for i := 4 to 8 do
writeln (10**i, ' samples give ', MC_Pi(i):7:5, ' as pi.');
end.
Output:
:> ./MonteCarlo 10000 samples give 3.14480 as pi. 100000 samples give 3.14484 as pi. 1000000 samples give 3.13970 as pi. 10000000 samples give 3.14100 as pi. 100000000 samples give 3.14162 as pi.
[edit] Perl
sub pi {
my $nthrows = shift;
my $inside = 0;
foreach (1 .. $nthrows) {
my $x = rand * 2 - 1,
$y = rand * 2 - 1;
if (sqrt($x*$x + $y*$y) < 1) {
$inside++;
}
}
return 4 * $inside / $nthrows;
}
printf "%9d: %07f\n", $_, pi($_) foreach 10**4, 10**6;
[edit] Perl 6
sub approximate_pi (Int $sample_size) {
my Int $in = 0;
(rand - 1/2)**2 + (rand - 1/2)**2 < 1/4 and ++$in
for ^$sample_size;
return 4 * $in / $sample_size;
}
say 'n = 100: ', approximate_pi 100;
say 'n = 1,000: ', approximate_pi 1_000;
say 'n = 10,000: ', approximate_pi 10_000;
And just to warp your brain, here's a version written in a purely functional style with no temp variables:
sub approximate_pi (Int $sample_size) {
$sample_size R/ [+]
4 xx grep 0 ..^ 1/4,
(rand - 1/2)**2 + (rand - 1/2)**2 xx
$sample_size;
}
The R metaoperator reverses the arguments to the division, producing the "goes into" operator.
Note also the use of a range for grep's test argument, which relies on the smartmatching semantics in Perl 6. Also, instead of counting by 1's and then multiplying by 4, we simply sum a series of 4's.
[edit] PHP
<?
$loop = 1000000; # loop to 1,000,000
$count = 0;
for ($i=0; $i<$loop; $i++) {
$x = rand() / getrandmax();
$y = rand() / getrandmax();
if(($x*$x) + ($y*$y)<=1) $count++;
}
echo "loop=".number_format($loop).", count=".number_format($count).", pi=".($count/$loop*4);
?>
Output:
loop=1,000,000, count=785,462, pi=3.141848
[edit] PicoLisp
(de carloPi (Scl)
(let (Dim (** 10 Scl) Dim2 (* Dim Dim) Pi 0)
(do (* 4 Dim)
(let (X (rand 0 Dim) Y (rand 0 Dim))
(when (>= Dim2 (+ (* X X) (* Y Y)))
(inc 'Pi) ) ) )
(format Pi Scl) ) )
(for N 6
(prinl (carloPi N)) )
Output:
3.4 3.23 3.137 3.1299 3.14360 3.140964
[edit] PowerShell
function Get-Pi ($Iterations = 10000) {
$InCircle = 0
for ($i = 0; $i -lt $Iterations; $i++) {
$x = Get-Random 1.0
$y = Get-Random 1.0
if ([Math]::Sqrt($x * $x + $y * $y) -le 1) {
$InCircle++
}
}
$Pi = [decimal] $InCircle / $Iterations * 4
$RealPi = [decimal] "3.141592653589793238462643383280"
$Diff = [Math]::Abs(($Pi - $RealPi) / $RealPi * 100)
New-Object PSObject `
| Add-Member -PassThru NoteProperty Iterations $Iterations `
| Add-Member -PassThru NoteProperty Pi $Pi `
| Add-Member -PassThru NoteProperty "% Difference" $Diff
}
This returns a custom object with appropriate properties which automatically enables a nice tabular display:
PS Home:\> 10,100,1e3,1e4,1e5,1e6 | ForEach-Object { Get-Pi $_ }
Iterations Pi % Difference
---------- -- ------------
10 3,6 14,591559026164641753596309630
100 3,40 8,225361302488828322840959090
1000 3,208 2,1138114877600474293158225800
10000 3,1444 0,0893606116311387583356211100
100000 3,14712 0,1759409006731298209938938800
1000000 3,141364 0,0072782698142600895432451100
[edit] PureBasic
OpenConsole()Output:
Procedure.d MonteCarloPi(throws.d)
inCircle.d = 0
For i = 1 To throws.d
randX.d = (Random(2147483647)/2147483647)*2-1
randY.d = (Random(2147483647)/2147483647)*2-1
dist.d = Sqr(randX.d*randX.d + randY.d*randY.d)
If dist.d < 1
inCircle = inCircle + 1
EndIf
Next i
pi.d = (4 * inCircle / throws.d)
ProcedureReturn pi.d
EndProcedure
PrintN ("'built-in' #Pi = " + StrD(#PI,20))
PrintN ("MonteCarloPi(10000) = " + StrD(MonteCarloPi(10000),20))
PrintN ("MonteCarloPi(100000) = " + StrD(MonteCarloPi(100000),20))
PrintN ("MonteCarloPi(1000000) = " + StrD(MonteCarloPi(1000000),20))
PrintN ("MonteCarloPi(10000000) = " + StrD(MonteCarloPi(10000000),20))
PrintN("Press any key"): Repeat: Until Inkey() <> ""
'built-in' #PI = 3.14159265358979310000 MonteCarloPi(10000) = 3.17119999999999980000 MonteCarloPi(100000) = 3.14395999999999990000 MonteCarloPi(1000000) = 3.14349599999999980000 MonteCarloPi(10000000) = 3.14127720000000020000 Press any key
[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
One use of the "sum" function is to count how many times something is true (because True = 1, False = 0):
>>> import random, math
>>> throws = 1000
>>> 4.0 * sum(math.hypot(*[random.random()*2-1
for q in [0,1]]) < 1
for p in xrange(throws)) / float(throws)
3.1520000000000001
>>> throws = 1000000
>>> 4.0 * sum(math.hypot(*[random.random()*2-1
for q in [0,1]]) < 1
for p in xrange(throws)) / float(throws)
3.1396359999999999
>>> throws = 100000000
>>> 4.0 * sum(math.hypot(*[random.random()*2-1
for q in [0,1]]) < 1
for p in xrange(throws)) / 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))
[edit] Faster implementation using Numpy
import numpy as np
n = input('Number of samples: ')
print np.sum(np.random.rand(n)**2+np.random.rand(n)<1)/float(n)*4
[edit] R
# nice but not suitable for big samples!
monteCarloPi <- function(samples) {
x <- runif(samples, -1, 1) # for big samples, you need a lot of memory!
y <- runif(samples, -1, 1)
l <- sqrt(x*x + y*y)
return(4*sum(l<=1)/samples)
}
# this second function changes the samples number to be
# multiple of group parameter (default 100).
monteCarlo2Pi <- function(samples, group=100) {
lim <- ceiling(samples/group)
olim <- lim
c <- 0
while(lim > 0) {
x <- runif(group, -1, 1)
y <- runif(group, -1, 1)
l <- sqrt(x*x + y*y)
c <- c + sum(l <= 1)
lim <- lim - 1
}
return(4*c/(olim*group))
}
print(monteCarloPi(1e4))
print(monteCarloPi(1e5))
print(monteCarlo2Pi(1e7))
[edit] Racket
#lang racket
(define (in-unit-circle? x y) (<= (sqrt (+ (sqr x) (sqr y))) 1))
;; point in ([-1,1], [-1,1])
(define (random-point-in-2x2-square) (values (* 2 (- (random) 1/2)) (* 2 (- (random) 1/2))))
;; Area of circle is (pi r^2). r is 1, area of circle is pi
;; Area of square is 2^2 = 4
;; There is a pi/4 chance of landing in circle
;; .: pi = 4*(proportion passed) = 4*(passed/samples)
(define (passed:samples->pi passed samples) (* 4 (/ passed samples)))
;; generic kind of monte-carlo simulation
(define (monte-carlo run-length report-frequency
sample-generator pass?
interpret-result)
(let inner ((samples 0) (passed 0) (cnt report-frequency))
(cond
[(= samples run-length) (interpret-result passed samples)]
[(zero? cnt) ; intermediate report
(printf "~a samples of ~a: ~a passed -> ~a~%"
samples run-length passed (interpret-result passed samples))
(inner samples passed report-frequency)]
[else
(inner (add1 samples)
(if (call-with-values sample-generator pass?)
(add1 passed) passed) (sub1 cnt))])))
;; (monte-carlo ...) gives an "exact" result... which will be a fraction.
;; to see how it looks as a decimal we can exact->inexact it
(let ((mc (monte-carlo 10000000 1000000 random-point-in-2x2-square in-unit-circle? passed:samples->pi)))
(printf "exact = ~a~%inexact = ~a~%(pi - guess) = ~a~%" mc (exact->inexact mc) (- pi mc)))
Output:
1000000 samples of 10000000: 785763 passed -> 785763/250000 2000000 samples of 10000000: 1571487 passed -> 1571487/500000 3000000 samples of 10000000: 2356776 passed -> 98199/31250 4000000 samples of 10000000: 3141924 passed -> 785481/250000 5000000 samples of 10000000: 3927540 passed -> 196377/62500 6000000 samples of 10000000: 4713072 passed -> 98189/31250 7000000 samples of 10000000: 5498300 passed -> 54983/17500 8000000 samples of 10000000: 6283199 passed -> 6283199/2000000 9000000 samples of 10000000: 7068065 passed -> 1413613/450000 exact = 3926793/1250000 inexact = 3.1414344 (pi - guess) = 0.00015825358979304482
A little more Racket-like is the use of an iterator (in this case for/fold), which is clearer than an inner function:
#lang racket
(define (in-unit-circle? x y) (<= (sqrt (+ (sqr x) (sqr y))) 1))
;; Good idea made in another task that:
;; The proportions of hits is the same in the unit square and 1/4 of a circle.
;; point in ([0,1], [0,1])
(define (random-point-in-unit-square) (values (random) (random)))
;; generic kind of monte-carlo simulation
;; Area of circle is (pi r^2). r is 1, area of circle is pi
;; Area of square is 2^2 = 4
;; There is a pi/4 chance of landing in circle
;; .: pi = 4*(proportion passed) = 4*(passed/samples)
(define (passed:samples->pi passed samples) (* 4 (/ passed samples)))
(define (monte-carlo/2 run-length report-frequency sample-generator pass? interpret-result)
(interpret-result
(for/fold ((pass 0))
([n (in-range run-length)]
#:when (when (and (not (zero? n)) (zero? (modulo n report-frequency)))
(printf "~a samples of ~a: ~a passed -> ~a~%"
n run-length pass (interpret-result pass n)))
#:when (call-with-values sample-generator pass?))
(add1 pass))
run-length))
;; (monte-carlo ...) gives an "exact" result... which will be a fraction.
;; to see how it looks as a decimal we can exact->inexact it
(let ((mc (monte-carlo/2 10000000 1000000 random-point-in-unit-square in-unit-circle? passed:samples->pi)))
(printf "exact = ~a~%inexact = ~a~%(pi - guess) = ~a~%" mc (exact->inexact mc) (- pi mc)))
[Similar output]
[edit] REXX
/*REXX program uses the Monte Carlo method to compute pi÷4 */
parse arg times chunks . /*does user want a specific num? */
if times=='' then times=1000000000 /*one billion should do it. */
if chunks=='' then chunks=10000 /*do Monte Carle in 10k chunks. */
limit=10000-1 /*REXX random gens only integers.*/
limitSq=limit**2 /*...so, instead of 1, use lim**2*/
!=0 /*number of "pi hits" so far. */
accur=0 /*accuracy of the Monte Carlo pi.*/
if 1=='f1'x then piChar='pi' /*if EBCDIC, then use literal. */
else piChar='e3'x /*if ASCII, then use symbol. */
pi=3.14159265358979323846264338327950288419716939937511 /*da real McCoy*/
numeric digits length(pi) /*at least, we'll use these digs.*/
say 'real pi='pi"+" /*might was well brag about it. */
say /*just for the eyeballs. */
do j=1 for times%chunks
do chunks /*do Monte Carlo, chunk-at-a-time*/
if random(0,limit)**2+random(0,limit)**2<=limitSq then !=!+1
end
reps=chunks*j /*compute number of repetitions. */
piX=4*!/reps /*let's see how this puppy does. */
_=compare(piX,pi) /*compare apples & ...crabapples.*/
if _<=accur then iterate /*if not better accuracy, pout. */
say right(comma(reps),20) 'repetitions: Monte Carlo' piChar,
"is accurate to" _-1 'places.' /*subtract 1 for dec point.*/
accur=_ /*use this accuracy for baseline.*/
end /*j*/
exit /*stick a fork in it, we're done.*/
/*────────────────────────────────COMMA subroutine──────────────────────*/
comma: procedure; parse arg _,c,p,t; arg ,cu; c=word(c ",",1)
if cu=='BLANK' then c=' '; o=word(p 3,1); p=abs(o); t=word(t 999999999,1)
if \datatype(p,'W') | \datatype(t,'W')|p==0|arg()>4 then return _; n=_'.9'
#=123456789; k=0; if o<0 then do; b=verify(_,' '); if b==0 then return _
e=length(_) - verify(reverse(_),' ') + 1; end; else do; b=verify(n,#,"M")
e=verify(n,#'0',,verify(n,#"0.",'M'))-p-1; end
do j=e to b by -p while k<t; _=insert(c,_,j); k=k+1; end; return _
output when using the default input:
real pi=3.14159265358979323846264338327950288419716939937511+
10,000 repetitions: Monte Carlo π is accurate to 3 places.
30,000 repetitions: Monte Carlo π is accurate to 4 places.
110,000 repetitions: Monte Carlo π is accurate to 5 places.
2,270,000 repetitions: Monte Carlo π is accurate to 6 places.
2,680,000 repetitions: Monte Carlo π is accurate to 7 places.
14,180,000 repetitions: Monte Carlo π is accurate to 9 places.
[edit] Ruby
def approx_pi(throws)
times_inside = throws.times.count {Math.hypot(rand, rand) <= 1.0}
# in pre-Ruby-1.8.7: times_inside = throws.times.select {Math.hypot(rand, rand) <= 1.0}.length
4.0 * times_inside / throws
end
[1000, 10_000, 100_000, 1_000_000, 10_000_000].each do |n|
puts "%8d samples: PI = %s" % [n, approx_pi(n)]
end
Output:
1000 samples: PI = 3.2 10000 samples: PI = 3.14 100000 samples: PI = 3.13244 1000000 samples: PI = 3.145124 10000000 samples: PI = 3.1414788
[edit] Scala
object MonteCarlo {
private val random = new scala.util.Random
/** Returns a random number between -1 and 1 */
def nextThrow: Double = (random.nextDouble * 2.0) - 1.0
/** Returns true if the argument point would be 'inside' the unit circle with
* center at the origin, and bounded by a square with side lengths of 2
* units. */
def insideCircle(pt: (Double, Double)): Boolean = pt match {
case (x, y) => (x * x) + (y * y) <= 1.0
}
/** Runs the simulation the specified number of times. Uses the result to
* estimate a value of pi */
def simulate(times: Int): Double = {
val inside = Iterator.tabulate (times) (_ => (nextThrow, nextThrow)) count insideCircle
inside.toDouble / times.toDouble * 4.0
}
def main(args: Array[String]): Unit = {
val sims = Seq(10000, 100000, 1000000, 10000000, 100000000)
sims.foreach { n =>
println(n+" simulations; pi estimation: "+ simulate(n))
}
}
}
Output:
10000 simulations; pi estimation: 3.1492 100000 simulations; pi estimation: 3.1396 1000000 simulations; pi estimation: 3.14208 10000000 simulations; pi estimation: 3.1409944 100000000 simulations; pi estimation: 3.1414386
[edit] Seed7
$ include "seed7_05.s7i";
include "float.s7i";
const func float: pi (in integer: throws) is func
result
var float: pi is 0.0;
local
var integer: throw is 0;
var integer: inside is 0;
begin
for throw range 1 to throws do
if rand(0.0, 1.0) ** 2 + rand(0.0, 1.0) ** 2 <= 1.0 then
incr(inside);
end if;
end for;
pi := flt(4 * inside) / flt(throws);
end func;
const proc: main is func
begin
writeln(" 10000: " <& pi( 10000) digits 5);
writeln(" 100000: " <& pi( 100000) digits 5);
writeln(" 1000000: " <& pi( 1000000) digits 5);
writeln(" 10000000: " <& pi( 10000000) digits 5);
writeln("100000000: " <& pi(100000000) digits 5);
end func;
Output:
10000: 3.14520
100000: 3.15000
1000000: 3.14058
10000000: 3.14223
100000000: 3.14159
[edit] Tcl
proc pi {samples} {
set i 0
set inside 0
while {[incr i] <= $samples} {
if {sqrt(rand()**2 + rand()**2) <= 1.0} {
incr inside
}
}
return [expr {4.0 * $inside / $samples}]
}
puts "PI is approx [expr {atan(1)*4}]\n"
foreach runs {1e2 1e4 1e6 1e8} {
puts "$runs => [pi $runs]"
}
result
PI is approx 3.141592653589793 1e2 => 2.92 1e4 => 3.1344 1e6 => 3.141924 1e8 => 3.14167724
[edit] Ursala
#import std
#import flo
mcp "n" = times/4. div\float"n" (rep"n" (fleq/.5+ sqrt+ plus+ ~~ sqr+ minus/.5+ rand)?/~& plus/1.) 0.
Here's a walk through.
-
mcp "n" =... defines a function namedmcpin terms of a dummy variable"n", which will be the number of iterations used in the simulation -
randignores its argument and returns a uniformly distributed number between 0 and 1 -
minus/.5is composed withrandto compute the difference between the random number and 0.5 -
sqrsquares the difference -
~~says to apply the function twice and return the pair of results -
pluscomposed with that adds the pair of results -
sqrttakes the square root of the sum -
fleq/.5is floating point comparison with a fixed right side of.5, returning true if its argument is greater or equal - Everything from
fleqtorandforms the predicate for the?conditional operator. - If the condition is true, the identity function is applied,
~& - If the condition is false, the
plus/1.function is applied, which adds one to its argument. -
rep"n"applied to a function has the effect of composing that function with itself"n"times, with"n"in this case being the parameter to themcpfunction. - The function being repeated
"n"times is applied to an argument of 0. - A division of the result by the number
"n"converted to a floating point value is performed bydiv\float"n". - The result of the division is quadrupled by
times/4..
test program:
#cast %eL
pis = mcp* <10,100,1000,10000,100000,1000000>
output:
< 2.800000e+00, 3.600000e+00, 3.164000e+00, 3.118800e+00, 3.144480e+00, 3.141668e+00>
[edit] XPL0
code Ran=1, CrLf=9;
code real RlOut=48;
func real MontePi(N); \Calculate pi using Monte Carlo method
int N; \number of randomly selected points
int I, X, Y, C;
def R = 10000; \radius of circle
[C:= 0; \initialize count of points in circle
for I:= 0 to N-1 do
[X:= Ran(R);
Y:= Ran(R);
if X*X + Y*Y <= R*R then C:= C+1;
];
return float(C)*4.0 / float(N); \Acir/Asqr = pi*R^2/4*R^2 = pi/4
];
[RlOut(0, MontePi( 100)); CrLf(0);
RlOut(0, MontePi( 10_000)); CrLf(0);
RlOut(0, MontePi( 1_000_000)); CrLf(0);
RlOut(0, MontePi(100_000_000)); CrLf(0);
]
Example output:
2.92000
3.13200
3.14375
3.14192
Categories:
- Programming Tasks
- Probability and statistics
- Ada
- ALGOL 68
- AutoHotkey
- AWK
- BASIC
- BBC BASIC
- C
- C sharp
- Clojure
- Common Lisp
- D
- E
- Euler Math Toolbox
- Factor
- Fantom
- Forth
- Fortran
- Go
- Haskell
- HicEst
- Icon
- Unicon
- Icon Programming Library
- J
- Java
- JavaScript
- K
- Liberty BASIC
- Locomotive Basic
- Logo
- LSL
- Lua
- Mathematica
- MATLAB
- MAXScript
- OCaml
- Octave
- PARI/GP
- Pascal
- Math
- Perl
- Perl 6
- PHP
- PicoLisp
- PowerShell
- PureBasic
- Python
- R
- Racket
- REXX
- Ruby
- Scala
- Seed7
- Tcl
- Ursala
- XPL0
