Monte Carlo methods
From Rosetta Code
(Redirected from Monte Carlo Simulation)
Monte Carlo methods is a programming task. Visitors like you are encouraged to solve it according to the task description, using any language they may happen to know.
A Monte Carlo Simulation is a way of approximating the value of a function where calculating the actual value is difficult or impossible. It uses random sampling to define constraints on the value and then makes a sort of "best guess."
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
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] ALGOL 68
Works with: ALGOL 68 version Standard - no extensions to language used
Works with: ALGOL 68G version Any - tested with release mk15-0.8b.fc9.i386
Works with: ELLA ALGOL 68 version Any (with appropriate job cards) - tested with release 1.8.8d.fc9.i386
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] 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] C
Translation of: Fortran
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
double drand(int lim)
{
/* there must be a better way (maybe) */
return ((double)lim * (double)rand() / (double)RAND_MAX );
}
double Pi(int samples)
{
int i, in_circle;
double coords[2], length;
in_circle = 0;
for(i=0; i<samples; i++)
{
coords[0] = drand(1);
coords[1] = drand(1);
coords[0] = coords[0]*2.0 - 1.0;
coords[1] = coords[1]*2.0 - 1.0;
length = sqrt(coords[0]*coords[0] + coords[1]*coords[1]);
if ( length <= 1.0 ) in_circle++;
}
return 4. * (double)in_circle / (double)samples;
}
int main()
{
int n = 10000;
while (n <= 100000000 )
{
printf("%d %lf\n", n, Pi(n));
n *= 10;
}
}
Output:
10000 3.119600 100000 3.143400 1000000 3.143064 10000000 3.142193 100000000 3.141625
[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: 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] 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] 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] 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 (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] 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);
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] 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] 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
[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] 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
Works with: Rakudo version #22 "Thousand Oaks"
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;
[edit] PowerShell
Works with: PowerShell version 2
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] 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] Ruby
def approx_pi(throws)
inside = throws.times.select {Math.hypot(rand, rand) <= 1.0}
4.0 * inside.length / 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] 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
Categories: Programming Tasks | Probability and statistics | Ada | ALGOL 68 | AutoHotkey | BASIC | C | C sharp | Clojure | Common Lisp | D | E | Factor | Forth | Fortran | Haskell | J | Java | Logo | Mathematica | MAXScript | OCaml | Octave | Perl | Perl 6 | PowerShell | PureBasic | Python | R | Ruby | Tcl
