% 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( library(random) ).
gamma(Mean,Alpha,Out) :-
Beta is Mean/Alpha,
(Alpha > 20
-> Mean is Alpha * Beta,
Sd is sqrt(Alpha*Beta*Beta),
normal(Mean,Sd,Out)
; gamma(Alpha,Beta,0,Out)).
gamma(0,_,X,X) :- !.
gamma(Alpha,Beta, In, Gamma) :-
random( R ),
Temp is In + ( -1 * Beta * log(1-R)),
Alpha1 is Alpha - 1,
gamma(Alpha1,Beta,Temp,Gamma).
% 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) :-
random( R1 ),
X1 is 2.0 * R1 - 1,
random( R2 ),
X2 is 2.0 * R2 - 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).