% 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).