:- module( hypergeometric, [hypergeometric/4] ).

:- ensure_loaded( library(lists) ).	% sum_list/2.
:- ensure_loaded( factorial ).		% /2.
:- pfd_demo:on_pl( sicstus(_), ensure_loaded(between) ).			% /3.
:- pfd_demo:on_pl( yap(_), ensure_loaded(between) ).			% /3.

% hypergeometric(  8, 32, 16 ).
% hypergeometric p80 hogg & tanis
% Collection of N = N1 + N2, objects of two distingushable classes
% I, (with N1 indistingushable members) and II (with N2 indistingushable
% objects). Choose R objects, without replacement. 
% hypergeometric/3 displays the distribution P(X=x), 0 >= x <= R, 
% which is the probability of x of the selected objects to have
% been of class I.
hypergeometric( N1, N2, R, Distro ) :-
	N is N1 + N2,
	bin_coef( N, R, Dnm ),
	findall( Px, 
		( between( 0, R, X ),
		  hypergeometric_point( X, R, N1, N2, Dnm, Px )
		),
		Distro ).

hypergeometric_point( X, R, N1, N2, Px ) :-
	N is N1 + N2,
	bin_coef( N, R, Dnm ),
	hypergeometric_point( X, R, N1, N2, Dnm, Px ).

hypergeometric_point( X, R, N1, N2, Dnm, Px ) :-
	RmX is R - X,
	bin_coef( N1, X, NmL ),
	bin_coef( N2, RmX, NmR ),
	IntNm is integer(NmL * NmR),
	IntDnm is integer(Dnm),
	Px  = IntNm / IntDnm.

% an older presentation 
hyper_old( N1, N2, R ) :-
	N is N1 + N2,
	bin_coef( N, R, Dnm ),
	hyper_old( 0, R, N1, N2, Dnm, 0, Sum ),
	write( sum:Sum ), nl.
hyper_old( X, R, N1, N2, Dnm, Acc, Sum ) :-
	( X > R -> 
		Sum = Acc 
		;
		RmX is R - X,
		bin_coef( N1, X, NmL ),
		bin_coef( N2, RmX, NmR ),
		Px is (NmL * NmR) / Dnm,
		write( X:Px ), nl,
		NxX is X + 1,
		NxAcc is Acc + Px,
		hyper_old( NxX, R, N1, N2, Dnm, NxAcc, Sum )
	).

% binomial_coefficients hogg & tanis, p 75.
bin_coef( N, R, Bcoef ) :-
	factorial( N, NFact ),
	factorial( R, RFact ),
	M is N - R,
	factorial( M, MFact ),
	Bcoef is NFact / RFact / MFact.