Monte Carlo methods

From Rosetta Code

(Redirected from Monte Carlo Simulation)
Jump to: navigation, search
Task
Monte Carlo methods
You are encouraged to solve this task according to the task description, using any language you may 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

D V.2. 

import std.stdio, std.random, std.math;
 
double pi(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 .. 10)
        writefln("%10s: %07f", 10^^p, pi(10^^p));
}

Sample output: 

        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] 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] Go

package main
import (
  "fmt"
  "rand"
  "math"
  "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.Nanoseconds())
  fmt.Println(getPi(10000))
  fmt.Println(getPi(100000))
  fmt.Println(getPi(1000000))
  fmt.Println(getPi(10000000))
  fmt.Println(getPi(100000000))
}

[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] 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] 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] MATLAB

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(abs(darts).^2,2) <= 1 );
 
    piEstimate = 4*sum(dartsInside)/numDarts;
 
end
 

Completely Vectorized: 

function piEstimate = monteCarloPi(numDarts)
 
    piEstimate = 4*sum(( sum(abs(rand(numDarts,2)).^2,2) <= 1 ))/numDarts;
 
end

Sample Output: 

>> monteCarloPi(30000000)
 
ans =
 
   3.141725200000000

[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] 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

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()
 
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() <> ""
 
Output: 
'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

[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 named mcp in terms of a dummy variable "n", which will be the number of iterations used in the simulation
  • rand ignores its argument and returns a uniformly distributed number between 0 and 1
  • minus/.5 is composed with rand to compute the difference between the random number and 0.5
  • sqr squares the difference
  • ~~ says to apply the function twice and return the pair of results
  • plus composed with that adds the pair of results
  • sqrt takes the square root of the sum
  • fleq/.5 is floating point comparison with a fixed right side of .5, returning true if its argument is greater or equal
  • Everything from fleq to rand forms 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 the mcp function.
  • 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 by div\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>
Retrieved from "http://rosettacode.org/wiki/Monte_Carlo_methods"
Personal tools