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