Monte Carlo methods: Difference between revisions

 

=={{header|11l}}==

V inside = 0

L 1..n

R 4.0 * inside / n

 

 

{{out}}

 

=={{header|360 Assembly}}==

MONTECAR CSECT

USING MONTECAR,R13 base register

WP DS PL16 packed decimal 16

YREGS

{{out}}

<pre>

{{libheader|Action! Tool Kit}}

{{libheader|Action! Real Math}}

 

DEFINE PTR="CARD"

Test(1000)

Test(10000)

{{out}}

[https://gitlab.com/amarok8bit/action-rosetta-code/-/raw/master/images/Monte_Carlo_methods.png Screenshot from Atari 8-bit computer]

 

=={{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>

=={{header|Arturo}}==

 

inside: new 0.0

do.times: throws [

loop [100 1000 10000 100000 1000000] 'n ->

 

{{out}}

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

=={{header|BASIC256}}==

{{works with|basic256|1.1.4.0}}

# Monte Carlo Simulator

# Determine value of pi

next i

 

{{out}}

<pre>

 

=={{header|BBC BASIC}}==

PRINT FNmontecarlo(10000)

PRINT FNmontecarlo(100000)

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

NEXT

{{out}}

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

 

=={{header|C++}}==

#include<iostream>

#include<math.h>

cout<<""<<4*double(hit)/double(imax)<<endl; } // Print out Pi number

}

 

=={{header|Clojure}}==

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

(if (< total iterations)

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

 

 

{{out}}

</pre>

 

[]

(if (<= (+ (Math/pow (rand) 2) (Math/pow (rand) 2)) 1) 1 0))

 

(pi-estimate 10000)

 

{{out}}

 

=={{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}}

=={{header|Crystal}}==

{{trans|Ruby}}

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

4.0 * times_inside / throws

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

puts "%8d samples: PI = %s" % [n, approx_pi(n)]

{{out}}

<pre> 1000 samples: PI = 3.1

 

=={{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

From example at [https://www.dartlang.org/ Dart Official Website]

 

import 'dart:async';

import 'dart:html';

bool get isInsideUnitCircle => x * x + y * y <= 1;

}

{{out}}

The script give in reality an output formatted in HTML

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:

 

=={{header|EasyLang}}==

for i range n

x = randomf

call mc 100000

call mc 1000000

Output:

3.1292

 

=={{header|Elixir}}==

def pi(n) do

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

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

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

 

{{out}}

=={{header|Erlang}}==

===With inline test===

-module(monte).

-export([main/1]).

 

main(N) -> io:format("PI: ~w~n", [ monte(N) ]).

{{out}}

<pre>

</pre>

===With test in a function===

-module(monte2).

-export([main/1]).

 

main(N) -> io:format("PI: ~w~n", [ monte(N) ]).

{{out}}

<pre>Xcoord

 

=={{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]]

There is some support and test expressions.

 

let print x = printfn "%A" x

 

MonteCarloPiGreco 10000 |> print

MonteCarloPiGreco 100000 |> print

{{out}}

<pre>

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

</pre>

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

program mc

integer :: n,i

pi = 4.d0 * dble( count( hypot(x(1,:),x(2,:)) <= 1.d0 ) ) / n

end function

 

=={{header|FreeBASIC}}==

' compile with: fbc -s console

 

Print : Print "hit any key to end program"

Sleep

{{out}}

<pre> Mumber of throws Ratio (Pi) Error

Since Futhark is a pure language, random numbers are implemented using Sobol sequences.

 

import "futlib/math"

 

let inside = reduce (+) 0 bs

in 4.0f32*f32(inside)/f32(n)

 

=={{header|Go}}==

'''Using standard library math/rand:'''

 

import (

fmt.Println(getPi(10000000))

fmt.Println(getPi(100000000))

{{out}}

<pre>

 

For very careful Monte Carlo studies, you might consider the subrepository rand library. The random number generator there has some advantages such as better known statistical properties and better use of memory.

 

import (

fmt.Println(getPi(10000000))

fmt.Println(getPi(100000000))

 

=={{header|Haskell}}==

import System.Random

 

let dist :: Double

dist = sqrt (rand_x * rand_x + rand_y * rand_y)

{{Out}}

<pre>Example:

Or, using foldM, and dropping sqrt:

 

import System.Random (randomRIO)

import Data.Functor ((<&>))

mapM_

(monteCarloPi >=> print)

{{Out}}

For example:

 

=={{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}}==

===ES5===

var x, y, m = 0;

 

console.log(mcpi(100000));

console.log(mcpi(1000000));

<pre>3.168

3.1396

 

===ES6===

"use strict";

 

);

});

{{Out}} For example:

<pre>1000 samples: 3.064

jq does not have a built-in PRNG so we will use /dev/urandom

as a source of entropy by invoking jq as follows:

function prng {

cat /dev/urandom | tr -cd '0-9' | fold -w 10 | sed 's/^0*\(.*\)*\(.\)*$/\1\2/'

}

 

 

'''program.jq'''

 

def percent: "\(100000 * . | round / 1000)%";

| ($p | mcPi) as $mcpi

| ((($pi - $mcpi)|length) / $pi) as $error

{{out}}

<pre>

=={{header|Jsish}}==

From Javascript ES5 entry, with PRNG seeded during unit testing for reproducibility.

function mcpi(n) {

var x, y, m = 0;

mcpi(1000000) ==> 3.142124

=!EXPECTEND!=

 

{{out}}

 

=={{header|Julia}}==

 

function monteπ(n)

p = monteπ(n)

println("$(lpad(n, 9)): π ≈ $(lpad(p, 10)), pct.err = ", @sprintf("%2.5f%%", abs(p - π) / π))

 

{{out}}

 

=={{header|K}}==

 

sim 10000

 

sim'10^!8

 

=={{header|Kotlin}}==

 

fun mcPi(n: Int): Double {

n *= 10

}

Sample output:

{{out}}

 

=={{header|Liberty BASIC}}==

for pow = 2 to 6

n = 10^pow

end function

 

{{out}}

=={{header|Locomotive Basic}}==

 

20 input "How many samples";n

30 u=n/100+1

140 print "Computed value of pi:"pi2!

150 print "Difference to real value of pi: ";

 

[[File:Monte Carlo, 200 points, Locomotive BASIC.png]]

 

=={{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

=={{header|Mathematica}}/{{header|Wolfram Language}}==

We define a function with variable sample size:

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

gives back:

<pre>10 3.2

10000000 3.14134</pre>

 

 

A less elegant way to solve the problem, is to imagine a (well-trained) monkey, throwing a number of darts at a dartboard.

 

We create a function ''MonkeyDartsPi'', which can take a variable number of throws as input:

xyCoordinates = RandomReal[{0, 1}, {numberOfThrows, 2}];

InsideCircle = Length[Select[Total[xyCoordinates^2, {2}],#<=1&]] ;

 

We do several runs with a larger number of throws each time, increasing by powers of 10.

 

We see that as the number of throws increases, we get closer to the value of Pi:

 

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 *

 

''Example:'' for n = ''1000'' the output is ''3.152''.

 

=={{header|Nim}}==

 

randomize()

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;

}

 

{{out}}

<pre>

 

=={{header|Phix}}==

<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>

<span style="color: #004080;">integer</span> <span style="color: #000000;">N</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">100</span>

<span style="color: #000000;">N</span> <span style="color: #0000FF;">*=</span> <span style="color: #000000;">10</span>

<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>

{{out}}

<pre>

 

=={{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>

Some general Monte Carlo simulators. <code>N</code> is the number of runs, <code>F</code> is the simulation function.

===Using while loop===

sim1(N, F) = C =>

C = 0,

C := C + apply(F),

I := I + 1

 

===List comprehension===

This is simpler, but slightly slower than using <code>while</code> loop.

 

===Recursion===

sim_rec(N,N,F,0,S).

sim_rec(0,_N,_F,S,S).

sim_rec(C,N,F,S0,S) :-

S1 = S0 + apply(F),

 

===Test===

Of the three different MC simulators, <code>sim_rec/2</code> (using recursion) is slightly faster than the other two (<code>sim1/2</code> and <code>sim2/2</code>) which have about the same speed.

foreach(N in 0..7)

sim_pi(10**N)

% The simulation function:

% returns 1 if success, 0 otherwise

 

{{out}}

 

=={{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 $_ }

 

=={{header|PureBasic}}==

Procedure.d MonteCarloPi(throws.d)

 

PrintN("Press any key"): Repeat: Until Inkey() <> ""

{{out}}

<pre>'built-in' #PI = 3.14159265358979310000

 

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|Quackery}}==

{{trans|Forth}}

 

 

[ [ 64 bit ] constant

swap 20 point$ echo$ cr ] is trials ( n --> )

 

 

{{out}}

 

=={{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]

{{works with|rakudo|2015-09-24}}

We'll consider the upper-right quarter of the unitary disk centered at the origin. Its area is <math>\pi \over 4</math>.

 

sub approximate_pi(Int $n) {

say "$_ iterations: ", approximate_pi $_

for 100, 1_000, 10_000;

{{out}}

<pre>Monte-Carlo π approximation:

We don't really need to write a function, though. A lazy list would do:

 

 

=={{header|REXX}}==

A specific─purpose commatizer function is included to format the number of iterations.

numeric digits 20 /*use 20 decimal digits to handle args.*/

parse arg times chunk digs r? . /*does user want a specific number? */

exit 0 /*stick a fork in it, we're all done. */

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

{{out|output|text=&nbsp; when using the default input:}}

<pre>

 

=={{header|Ring}}==

decimals(8)

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

t = (4 * n) / t

return t

Output:

<pre>

 

=={{header|Ruby}}==

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

4.0 * times_inside / throws

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

puts "%8d samples: PI = %s" % [n, approx_pi(n)]

{{out}}

<pre> 1000 samples: PI = 3.2

 

=={{header|Rust}}==

 

use rand::Rng;

println!("{:9}: {:<11} dev: {:.5}%", samples, estimate, deviation);

}

{{out}}

<pre>Real pi: 3.141592653589793

 

=={{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}}

=={{header|SequenceL}}==

First solution is serial due to the use of random numbers. Will always give the same result for a given n and seed

import <Utilities/Random.sl>;

import <Utilities/Conversion.sl>;

result when n < 0 else

monteCarloHelper(n - 1, yRand.Generator, newResult);

 

The second solution will run in parallel. It will also always give the same result for a given n and seed. (Note, the function monteCarloHelper is the same in both versions).

 

import <Utilities/Random.sl>;

import <Utilities/Conversion.sl>;

result when n < 0 else

monteCarloHelper(n - 1, yRand.Generator, newResult);

 

=={{header|Sidef}}==

4 * (^nthrows -> count_by {

hypot(1.rand(2) - 1, 1.rand(2) - 1) < 1

for n in [1e2, 1e3, 1e4, 1e5, 1e6] {

printf("%9d: %07f\n", n, monteCarloPi(n))

{{out}}

<pre>

 

=={{header|Stata}}==

clear all

quietly set obs `1'

 

. mcdisk 100000000

 

=={{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><

{{trans|Kotlin}}

{{libheader|Wren-fmt}}

import "/fmt" for Fmt

 

Fmt.print("$9d -> $10.8f -> $6.4f", n, pi, err)

n = n * 10

 

{{out}}

 

=={{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>