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