:- 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.