Monte Carlo methods: Difference between revisions

{{task|Probability and statistics}}

A '''Monte Carlo Simulation''' is a way of approximating the value of a function

where calculating the actual value is difficult or impossible. <br>

and then makes a sort of "best guess."

 

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

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

 

Also, show the results of a few different sample sizes.

 

 

=={{header|Ada}}==

with Ada.Numerics.Float_Random; use Ada.Numerics.Float_Random;

 

Put_Line (" 10_000_000:" & Float'Image (Pi ( 10_000_000)));

Put_Line ("100_000_000:" & Float'Image (Pi (100_000_000)));

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.

 

{{works with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release 1.8.8d.fc9.i386}}

BEGIN

INT inside := 0;

print ((" 1 000 000:",pi ( 1 000 000),new line));

print ((" 10 000 000:",pi ( 10 000 000),new line));

{{out}}

<pre>

</pre>

 

{{AutoHotkey case}}

Source: [http://www.autohotkey.com/forum/topic44657.html AutoHotkey forum] by Laszlo

MsgBox % MontePi(10000) ; 3.154400

MsgBox % MontePi(100000) ; 3.142040

Return 4*p/n

}

 

=={{header|AWK}}==

# --- with command line argument "throws" ---

 

Pi = 3.14333

 

 

=={{header|BASIC}}==

{{works with|QuickBasic|4.5}}

{{trans|Java}}

CLS

PRINT getPi(10000)

NEXT i

getPi = 4! * inCircle / throws

{{out}}

<pre>

</pre>

 

PRINT FNmontecarlo(10000)

PRINT FNmontecarlo(100000)

IF RND(1)^2 + RND(1)^2 < 1 n% += 1

NEXT

{{out}}

<pre>

3.1412816

</pre>

 

=={{header|C}}==

#include <stdlib.h>

#include <math.h>

printf("Pi is %f\n", pi(3e-4)); /* set to 1e-4 for some fun */

return 0;

 

=={{header|C sharp|C#}}==

 

class Program {

}

}

 

{{out}}

100,000,000:3.1413976

</pre>

 

=={{header|Clojure}}==

(loop [x (rand) y (rand) in 0 total 1]

(if (< total iterations)

(double (* (/ in total) 4)))))

 

 

{{out}}

100000: 3.14104

1000000: 3.141064

</pre>

 

=={{header|Common Lisp}}==

(/ (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))

 

{{out}}

1000000 -> 3.142072

10000000 -> 3.1420677

</pre>

 

=={{header|D}}==

 

double pi(in uint nthrows) /*nothrow*/ @safe /*@nogc*/ {

foreach (immutable p; 1 .. 8)

writefln("%10s: %07f", 10 ^^ p, pi(10 ^^ p));

{{out}}

<pre> 10: 3.200000

 

===More Functional Style===

import std.stdio, std.random, std.math, std.algorithm, std.range;

 

foreach (immutable p; 1 .. 8)

writefln("%10s: %07f", 10 ^^ p, pi(10 ^^ p));

{{out}}

<pre> 10: 3.200000

1000000: 3.142836

10000000: 3.141550</pre>

 

=={{header|E}}==

This computes a single quadrant of the described square and circle; the effect should be the same since the other three are symmetric.

 

var inside := 0

for _ ? (entropy.nextFloat() ** 2 + entropy.nextFloat() ** 2 < 1) in 1..n {

}

return inside * 4 / n

 

Some sample runs:

? pi(100000)

# value: 3.13848

 

=={{header|Elixir}}==

def pi(n) do

count = Enum.count(1..n, fn _ ->

:math.sqrt(x*x + y*y) <= 1

end)

Enum.each([1000, 10000, 100000, 1000000, 10000000], fn n ->

:io.format "~8w samples: PI = ~f~n", [n, MonteCarlo.pi(n)]

 

{{out}}

1000000 samples: PI = 3.142904

10000000 samples: PI = 3.141124

</pre>

 

=={{header|ERRE}}==

PROGRAM RANDOM_PI

 

MONTECARLO(1000000->RES) PRINT(RES)

MONTECARLO(10000000->RES) PRINT(RES)

{{out}}

<pre>

 

=={{header|Euler Math Toolbox}}==

>function map MonteCarloPI (n,plot=false) ...

$ X:=random(1,n);

3.14159265359

>MonteCarloPI(10000,true):

 

[[File:Test.png]]

 

=={{header|Factor}}==

Since Factor lets the user choose the range of the random generator, we use 2^32.

 

 

: limit ( -- n ) 2 32 ^ ; inline

: in-circle ( x y -- ? ) limit [ sq ] tri@ [ + ] [ <= ] bi* ;

: rand ( -- r ) limit random ;

 

Example use:

 

 

=={{header|Fantom}}==

 

class MontyCarlo

{

}

}

 

{{out}}

=={{header|Fortran}}==

{{works with|Fortran|90 and later}}

IMPLICIT NONE

END DO

 

{{out}}

100000000 3.14147

</pre>

 

=={{header|Go}}==

 

import (

fmt.Println(getPi(10000000))

fmt.Println(getPi(100000000))

{{out}}

<pre>

3.14149596

</pre>

 

=={{header|Haskell}}==

import System.Random

 

where

Prelude System.Random Control.Monad> get_pi 10000

3.1352

3.15184

Prelude System.Random Control.Monad> get_pi 1000000

 

=={{header|HicEst}}==

inside = 0

DO i = 1, samples

WRITE(ClipBoard) Pi(1E5) ! 3.14204

WRITE(ClipBoard) Pi(1E6) ! 3.141672

 

=={{header|Icon}} and {{header|Unicon}}==

every t := 10 ^ ( 5 to 9 ) do

printf("Rounds=%d Pi ~ %r\n",t,getPi(t))

incircle +:= 1

return 4 * incircle / rounds

 

{{libheader|Icon Programming Library}}

=={{header|J}}==

'''Explicit Solution:'''

4* y%~ +/ 1>: %: +/ *: <: +: (2,y) ?@$ 0

 

'''Tacit Solution:'''

 

'''Examples:'''

3.1426

piMC 10^i.7

 

=={{header|Java}}==

public static void main(String[] args) {

System.out.println(getPi(10000));

return 4.0 * inCircle / numThrows;

}

{{out}}

3.1396

3.14168604

{{works with|Java|8+}}

 

import java.util.stream.IntStream;

return (4.0 * inCircle) / numThrows;

}

{{out}}

3.1556

 

=={{header|JavaScript}}==

 

 

 

}

 

console.log(mcpi(100000));

console.log(mcpi(1000000));

<pre>3.168

3.1396

</pre>

 

 

{{out}}

<pre>

</pre>

 

=={{header|K}}==

 

sim 10000

 

sim'10^!8

 

 

 

{{out}}

<pre>

</pre>

 

=={{header|Logo}}==

to square :n

output :n * :n

show sim 100000 10000 ; 3.145

show sim 1000000 10000 ; 3.140828

 

=={{header|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 iMAX_SAMPLE_POWER = 6;

default {

llOwnerSay("Done.");

}

{{out}}

<pre>Estimating Pi (3.141593)

 

=={{header|Lua}}==

math.randomseed( os.time() )

 

print( MonteCarlo( 100000 ) )

print( MonteCarlo( 1000000 ) )

{{out}}

<pre>3.1436

3.1420188</pre>

 

We define a function with variable sample size:

Example (samplesize=10,100,1000,....10000000):

gives back:

100 3.16

1000 3.152

100000 3.14872

1000000 3.1408

 

 

=={{header|MATLAB}}==

 

Minimally Vectorized:

 

%The square has a sides of length 2, which means the circle has radius

 

end

 

Completely Vectorized:

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

 

 

{{out}}

 

ans =

 

=={{header|Maxima}}==

approx_pi(n):= block(

[x: random_continuous_uniform(0, 1, n),

4*cin/n);

 

=={{header|MAXScript}}==

 

=={{header|МК-61/52}}==

1 - x<0 15 КИП4 L0 04 ИП4 4 *

 

 

=={{header|Nim}}==

 

for i in 1..int64(nthrows):

for n in [10e4, 10e6, 10e7, 10e8]:

{{out}}

<pre>3.15336

 

=={{header|OCaml}}==

let rec helper i count =

if i = throws then count

else

helper (i+1) count

Example:

# get_pi 10000;;

=={{header|Octave}}==

 

in_circle = 0;

for samp = 1:samples

disp(montepi(l));

l *= 10;

 

Since it runs slow, I've stopped it at the second iteration, obtaining:

=== Much faster implementation ===

 

function result = montepi(n)

result = sum(rand(1,n).^2+rand(1,n).^2<1)/n*4;

endfunction

 

=={{header|PARI/GP}}==

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.

 

=={{header|Pascal}}==

{{libheader|Math}}

 

uses

writeln (10**i, ' samples give ', MC_Pi(i):7:5, ' as pi.');

end.

{{out}}

<pre>:> ./MonteCarlo

 

=={{header|Perl}}==

my $nthrows = shift;

my $inside = 0;

foreach (1 .. $nthrows) {

if (sqrt($x*$x + $y*$y) < 1) {

$inside++;

}

 

 

{{out}}

 

=={{header|PHP}}==

$loop = 1000000; # loop to 1,000,000

$count = 0;

}

echo "loop=".number_format($loop).", count=".number_format($count).", pi=".($count/$loop*4);

{{out}}

<pre>loop=1,000,000, count=785,462, pi=3.141848</pre>

 

=={{header|PicoLisp}}==

(let (Dim (** 10 Scl) Dim2 (* Dim Dim) Pi 0)

(do (* 4 Dim)

 

(for N 6

{{out}}

<pre>3.4

=={{header|PowerShell}}==

{{works with|PowerShell|2}}

$InCircle = 0

for ($i = 0; $i -lt $Iterations; $i++) {

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

<pre>PS Home:\> 10,100,1e3,1e4,1e5,1e6 | ForEach-Object { Get-Pi $_ }

100000 3,14712 0,1759409006731298209938938800

1000000 3,141364 0,0072782698142600895432451100</pre>

 

=={{header|Python}}==

 

One use of the "sum" function is to count how many times something is true (because True = 1, False = 0):

>>> throws = 1000

>>> 4.0 * sum(math.hypot(*[random.random()*2-1

for q in [0,1]]) < 1

for p in xrange(throws)) / float(throws)

 

===As a program using a function===

from random import random

from math import hypot

for n in [10**4, 10**6, 10**7, 10**8]:

print "%9d: %07f" % (n, pi(n))

 

===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)**2<1)/float(n)*4

 

=={{header|R}}==

monteCarloPi <- function(samples) {

x <- runif(samples, -1, 1) # for big samples, you need a lot of memory!

print(monteCarloPi(1e4))

print(monteCarloPi(1e5))

 

=={{header|Racket}}==

 

(define (in-unit-circle? x y) (<= (sqrt (+ (sqr x) (sqr y))) 1))

;; 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)))

{{out}}

<pre>1000000 samples of 10000000: 785763 passed -> 785763/250000

A little more Racket-like is the use of an iterator (in this case '''for/fold'''),

which is clearer than an inner function:

(define (in-unit-circle? x y) (<= (sqrt (+ (sqr x) (sqr y))) 1))

;; Good idea made in another task that:

;; 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)))

 

[Similar output]

 

=={{header|REXX}}==

/*──────────────────────────────────────────────────────────────────────────────────────*/

<pre>

 

 

</pre>

 

=={{header|Ring}}==

decimals(8)

see "monteCarlo(1000) = " + monteCarlo(1000) + nl

t = (4 * n) / t

return t

Output:

<pre>

monteCarlo(10000) = 3.00320000

monteCarlo(100000) = 2.99536000

</pre>

 

=={{header|Ruby}}==

times_inside = throws.times.count {Math.hypot(rand, rand) <= 1.0}

4.0 * times_inside / throws

end

 

[1000, 10_000, 100_000, 1_000_000, 10_000_000].each do |n|

{{out}}

<pre> 1000 samples: PI = 3.2

1000000 samples: PI = 3.145124

10000000 samples: PI = 3.1414788</pre>

 

=={{header|Scala}}==

private val random = new scala.util.Random

 

}

}

{{out}}

<pre>10000 simulations; pi estimation: 3.1492

 

=={{header|Seed7}}==

include "float.s7i";

 

writeln(" 10000000: " <& pi( 10000000) digits 5);

writeln("100000000: " <& pi(100000000) digits 5);

 

{{out}}

100000000: 3.14159

</pre>

=={{header|Swift}}==

{{trans|JavaScript}}

 

func mcpi(sampleSize size:Int) -> Double {

println(mcpi(sampleSize: 1000000))

println(mcpi(sampleSize: 10000000))

{{out}}

<pre>

 

=={{header|Tcl}}==

set i 0

set inside 0

foreach runs {1e2 1e4 1e6 1e8} {

puts "$runs => [pi $runs]"

result

<pre>PI is approx 3.141592653589793

 

=={{header|Ursala}}==

#import flo

 

Here's a walk through.

* <code>mcp "n" = </code>... defines a function named <code>mcp</code> in terms of a dummy variable <code>"n"</code>, which will be the number of iterations used in the simulation

* The result of the division is quadrupled by <code>times/4.</code>.

test program:

 

{{out}}

<pre><

3.144480e+00,

3.141668e+00></pre>

 

=={{header|XPL0}}==

code real RlOut=48;

 

RlOut(0, MontePi( 1_000_000)); CrLf(0);

RlOut(0, MontePi(100_000_000)); CrLf(0);

 

{{out}}

 

=={{header|zkl}}==

inCircle:=0;

do(n){

}

4.0*inCircle/n

Or, in a more functional style (using a reference for state info):

4.0 * (1).pump(n,Void,fcn(r){

x:=(0.0).random(1); y:=(0.0).random(1);

r

}.fp(Ref(0)) ).value/n;

{{out}}

<pre>