|
|
Random numbers in Prolog
Rand-Num: generating random variables in Prolog
ranf.pl=
% ranf.pl
% random floats 0 <= x < 1
% the swi-prolog "is" function can
% be programmed. predicates of
% arity N can be called as functions
% of arity N-1
:- arithmetic_function(ranf/0).
ranf(X) :-
N = 65536,
% X is random(N) returns a
% random number from 0 .. N-1
X is random(N)/N.
|
normal.pl=
% normal.pl
% using ranf to generate a normal distribution
% -inf <= x <= inf, mean=M, standard deviation=S
% if "ranf" already loaded,
% then nothing happens
:- ensure_loaded(ranf).
:- arithmetic_function(normal/2).
normal(M,S,N) :-
box_muller(M,S,N).
% classical fast method for computing
% normal functions using polar co-ords
% (no- i dont understand it either)
box_muller(M,S,N) :-
w(W0,X),
W is sqrt((-2.0 * log(W0))/W0),
Y1 is X * W,
N is M + Y1*S.
w(W,X) :-
X1 is 2.0 * ranf - 1,
X2 is 2.0 * ranf - 1,
W0 is X1*X1 + X2*X2,
% IF -> THEN ; ELSE
% same as xx :- IF,!, THEN,
% xx :- ELSE
% vars bound in IF not available to ELSE
% no backtracking within the IF
% -> ; precendence higher than ,
(W0 >= 1.0 -> w(W,X) ; W0=W, X = X1).
|
beta.pl=
% beta.pl
% simple beta distributions (0 <= x <=1, mean=B)
% only B = 0.1, 0.2, 0.3, ... 0.9, 1 is supported
:- ensure_loaded(ranf).
:- arithmetic_function(beta/1).
beta(B,X) :- beta1(B,X),!.
beta(B,X) :- B1 is 1 - B, beta1(B1,Y),X is 1 - Y.
% the following tricks came from a statistics
% guy who did the algebra and showed me some
% tricks for quickly generating beta vars using
% fractional powers.
beta1(0.50,X) :- X is ranf.
beta1(0.60,X) :- X is ranf^0.67.
beta1(0.67,X) :- X is ranf^0.5.
beta1(0.75,X) :- X is ranf^0.33.
beta1(0.80,X) :- X is ranf^0.25.
beta1(0.9,X) :- X is ranf^(1/9).
beta1(1,1). |
mean = a B
variance = a B2
gamma.pl=
% gamma.pl
% 0 <= x <= inf
% mean = Mean, skew= Alpha
% at low skew, (e.g. x=2), x always near mean
% at skew=20, gamma becomes normal
% only supports integer values X and Alpha
:- ensure_loaded(normal).
:- arithmetic_function(gamma/2).
gamma(Mean,Alpha,Out) :-
Beta is Mean/Alpha,
(Alpha > 20
-> Mean is Alpha * Beta,
Sd is sqrt(Alpha*Beta*Beta),
Out is normal(Mean,Sd,Out)
; gamma(Alpha,Beta,0,Out)).
gamma(0,_,X,X) :- !.
gamma(Alpha,Beta, In, Gamma) :-
Temp is In + ( -1 * Beta * log(1-ranf)),
Alpha1 is Alpha - 1,
gamma(Alpha1,Beta,Temp,Gamma). |
random.pl=
% random.pl
:- ensure_loaded(format),
ensure_loaded(bins),
ensure_loaded(ranf),
ensure_loaded(normal),
ensure_loaded(beta),
ensure_loaded(gamma).
demo :- tell('random.out'), ignore(demo1), told.
demo1 :-
Repeats is 10^4,
forall(member(F, [normal(20,2),
10*beta(0.8),
gamma(10,15),
gamma(10,5),
gamma(10,2)
]),
demo2(Repeats,F)).
demo2(Repeats,F) :-
format('\n\n---| ~w * ~w |-------',
[Repeats,F]),
Size=1,
% standard "REPEATS times do"
findall(X,(between(1,Repeats,_)
% here's a wierd one: calling
% the function passed in as data
,X is F
),L0),
cutDown2Sizes(Size,L0,L),
dist(5,5,100,L).
cutDown2Sizes(Size) -->
maplist(cutDown2Size(Size)).
cutDown2Size(Size,X,Y) :-
Y is round(X/Size).
|
Which generates
random.out=
---| 10000 * normal(20, 2) |-------
item frequency
28 2
27 6
26 23
25 100 ~
24 287 ~~~
23 636 ~~~~~~
22 1240 ~~~~~~~~~~~~
21 1730 ~~~~~~~~~~~~~~~~~
20 1944 ~~~~~~~~~~~~~~~~~~~
19 1754 ~~~~~~~~~~~~~~~~~~
18 1243 ~~~~~~~~~~~~
17 628 ~~~~~~
16 274 ~~~
15 98 ~
14 28
13 7
---| 10000 * 10*beta(0.8) |-------
item frequency
10 1841 ~~~~~~~~~~~~~~~~~~
9 2926 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
8 2060 ~~~~~~~~~~~~~~~~~~~~~
7 1353 ~~~~~~~~~~~~~~
6 896 ~~~~~~~~~
5 503 ~~~~~
4 271 ~~~
3 110 ~
2 36
1 4
---| 10000 * gamma(10, 15) |-------
item frequency
22 1
21 1
20 7
19 15
18 43
17 96 ~
16 141 ~
15 239 ~~
14 425 ~~~~
13 673 ~~~~~~~
12 983 ~~~~~~~~~~
11 1257 ~~~~~~~~~~~~~
10 1482 ~~~~~~~~~~~~~~~
9 1550 ~~~~~~~~~~~~~~~~
8 1391 ~~~~~~~~~~~~~~
7 971 ~~~~~~~~~~
6 507 ~~~~~
5 171 ~~
4 45
3 2
---| 10000 * gamma(10, 5) |-------
item frequency
38 1
32 2
31 2
30 1
29 1
28 7
27 7
26 16
25 19
24 24
23 40
22 53 ~
21 82 ~
20 105 ~
19 126 ~
18 174 ~~
17 201 ~~
16 294 ~~~
15 383 ~~~~
14 461 ~~~~~
13 586 ~~~~~~
12 678 ~~~~~~~
11 735 ~~~~~~~
10 867 ~~~~~~~~~
9 951 ~~~~~~~~~~
8 965 ~~~~~~~~~~
7 948 ~~~~~~~~~
6 793 ~~~~~~~~
5 682 ~~~~~~~
4 449 ~~~~
3 251 ~~~
2 87 ~
1 9
---| 10000 * gamma(10, 2) |-------
item frequency
56 1
52 1
51 1
49 4
48 3
47 1
46 3
45 3
44 6
43 5
42 3
41 6
40 4
39 5
38 2
37 11
36 6
35 9
34 19
33 12
32 20
31 28
30 23
29 36
28 35
27 58 ~
26 57 ~
25 67 ~
24 79 ~
23 96 ~
22 107 ~
21 134 ~
20 162 ~~
19 162 ~~
18 221 ~~
17 236 ~~
16 276 ~~~
15 284 ~~~
14 332 ~~~
13 377 ~~~~
12 429 ~~~~
11 479 ~~~~~
10 546 ~~~~~
9 628 ~~~~~~
8 607 ~~~~~~
7 684 ~~~~~~~
6 666 ~~~~~~~
5 765 ~~~~~~~~
4 708 ~~~~~~~
3 668 ~~~~~~~
2 552 ~~~~~~
1 328 ~~~
0 45
|
Support code for the above
bins.pl=
% bins.pl
% code to generate histograms. cool just for
% its succinctness
:- ensure_loaded(format).
bins(L0,L) :-
% "sort" removes duplicates
% "msort" keeps copies
msort(L0,L1),
bins(L1,[],L).
% returns a list showing how often a
% particular item appears. assumes the
% input list is sorted
bins([],X,X).
bins([H|T],[H-N0|Rest],Out) :- !,
N is N0 + 1,
bins(T,[H-N|Rest],Out).
bins([H|T],In,Out) :-
bins(T,[H-1|In],Out).
dist(L) :- dist(5,5,3,L).
dist(W1, % width of the first "item" column
W2, % width of the second "frequency" column
W3, % scale factor for how much to shrink
% the twiddles shown in column three
L) :-
% sformat builds a string and binds it to "S".
% this string stores the widths and scale factot
% for our columns. note the use of ">" and "S"
% which are special format commands defined in
% format.pl
sformat(S,'~~~w> ~~~w> ~~~wS\n',[W1,W2,W3]),
bins(L,Bins),
nl,
format(S,[item,frequency,0]),
% so succint: just loop through the inputs, calling
% format on each item using our special "S" string
forall(member(What-N,Bins),format(S,[What,N,N])). |
format.pl=
% format.pl
%--------------------------------------
% stuff to simplify printing clauses.
% swi-prolog lets a programmer customize
% the format statement.
% FIRST, a predicate of arity 2 is registerred
% next to some letter
:- format_predicate('P',p(_,_)).
% SECOND, write the predicate.
p(default,X) :- !, p(0,X).
p(_,(X :- true)) :- !, format('~p.\n',X).
p(_,(X :- Y )) :- !, portray_clause((X :- Y)).
p(N,[H|T] ) :- !, not(not((numbervars([H|T],N,_),
format('~p',[[H|T]])))).
p(N,X ) :- not(not((numbervars(X,N,_),
format('~p',X)))).
%--------------------------------------
% stuff to simplify right justifying text
% FIRST:
:- format_predicate('>',padChars(_,_)).
% SECOND:
padChars(default,A) :-
padChars(5,A).
% the first arg "S" is the optional argument
% someone may have given with the "~" command
padChars(S,A) :-
writeThing(A,Thing,N),
Pad is S - N,
% standard trick to emulate
% for(i=1;i<=N;i++) { doThis }
forall(between(1,Pad,_),put(32)),
write(Thing).
writeThing(X,S,L) :-
% sformat returns the string in
% the first arg
sformat(S,'~w',[X]), string_length(S,L).
%--------------------------------------
% stuff to simplify left justifying text
% FIRST:
:- format_predicate('<',charsPad(_,_)).
% SECOND:
charsPad(default,A) :- charsPad(5,A).
charsPad(S,A) :-
writeThing(A,Thing,N),
atom_length(A,N),
Pad is S - N,
write(Thing),
forall(between(1,Pad,_),put(32)).
%--------------------------------------
% stuff to simplify printing N twiddles,
% scaled to some factor
% FIRST:
:- format_predicate('S',twiddle(_,_)).
% SECOND:
twiddle(default,A) :- twiddle(25,A).
twiddle(W,N) :-
N1 is round(N/W),
forall(between(1,N1,_),put(126)).
%--------------------------------------
% stuff to simplify printing lists
% scaled to some factor
% FIRST:
:- format_predicate('L',printL(_,_)).
% SECOND:
printL(default,List) :- printL(10^6,List).
printL(TooLong,List) :-
forall((nth1(Pos,List,Item),
Pos < TooLong),
format('\t~w\n',Item)).
|
Not © Tim Menzies, 2001
Share and enjoy- information wants to be free.
But if you take anything from this site,
please credit tim@menzies.com.
|
|
