/
monte_carlo_methods.pi
84 lines (64 loc) · 2.18 KB
/
monte_carlo_methods.pi
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
/*
Monte Carlo methods in Picat.
From Rosetta code:
http://rosettacode.org/wiki/Monte_Carlo_methods
"""
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
"""
Model created by Hakan Kjellerstrand, hakank@gmail.com
See also my Picat page: http://www.hakank.org/picat/
*/
import math.
main => go.
go =>
foreach(N in 0..6)
time(sim_pi(10**N))
end,
nl.
go2 =>
foreach(N in 0..6)
time(sim_pi2(10**N))
end,
nl.
sim_pi(N) =>
println($sim_pi(N)),
Inside = sim(N, pi_f),
println(inside=Inside),
MyPi = 4.0*Inside/N,
println(myPi=MyPi),
Pi = 4*atan2(1,1),
println(pi=Pi),
println([n=N, myPi=MyPi, diff=Pi-MyPi]).
sim_pi2(N) =>
Inside = sim2(N, pi_f),
MyPi = 4.0*Inside/N,
Pi = 4*atan2(1,1),
println([n=N, myPi=MyPi, diff=Pi-MyPi]).
% the simulation function:
% returns 1 if success, 0 otherwise
pi_f() = cond(frand()**2 + frand()**2 <= 1, 1, 0).
% simple (but general) Monte Carlo simulator
% f() is the simulation function
sim(N, F) = C =>
println($sim(N, F)),
C = 0,
foreach (_I in 1..N)
C := C + apply($F)
end.
sim2(N, F) = sum([apply(F) : _I in 1..N]).
% random value between 0.0..1.0
randf() = V => (_I,V) = modf(random()/1000000.0).