1:- module(prob_tagged,
    2          [ uniform01//1
    3          , uniform//2, uniformT//2, uniformP//2
    4          , normal//1
    5          , gaussian//3
    6          , exponential//1
    7          , poisson//2
    8          , stable//3
    9          , dirichlet//2
   10          , discrete//2
   11          , discrete//3, discreteT//3
   12          , binomial//3
   13          , beta//3
   14          , zeta//2
   15          , gamma//2
   16          , inv_gamma//2
   17          , bernoulli//2
   18          , students_t//2
   19          , mixture//3
   20          , product_pair//3
   21          , prob/3, prob/2
   22          , pdf/3, pdf/2
   23          ]).

Random predicates

This module provides a set of predicates for sampling from various distributions. The state of the random generator is threaded through using the DCG idiom.

author
- Samer Abdallah /
   34:- module_transparent stream/2.   35
   36:- use_module(library(dcg_core)).   37:- use_module(library(dcg_pair)).   38:- use_module(library(plrand),[]).   39
   40term_expansion(wrap_rs(Arity,Name,Pred), Head :- plrand:Body) :-
   41   length(A1, Arity),
   42   append(A1, [rs(S1),rs(S2)], HeadArgs), Head =.. [Name | HeadArgs],
   43   append(A1, [S1,S2],         PredArgs), Body =.. [Pred | PredArgs].
   44
   45term_expansion(wrap_disc(Arity,Name,Pred), Head :- (plrand:Body, P2 is P1*P)) :-
   46   length(A1, Arity),
   47   append(A1, [p(P1),p(P2)], HeadArgs), Head =.. [Name | HeadArgs],
   48   append(A1, [P],           PredArgs), Body =.. [Pred | PredArgs].
   49
   50term_expansion(wrap_cont(Arity,Name,Pred), [Head1, (Head :- (plrand:Body, P2 is P1*P))]) :-
   51   length(A1, Arity),
   52   append(A1, [p(_),p(0)],     Head1Args), Head1 =.. [Name | Head1Args],
   53   append(A1, [pd(P1),pd(P2)], HeadArgs), Head =.. [Name | HeadArgs],
   54   append(A1, [P],             PredArgs), Body =.. [Pred | PredArgs].
 bernoulli(+A:prob, -X:oneof([0,1]))// is det
Sample binary random variable.
   58bernoulli(P,X,rs(S1),rs(S2)) :- !, plrand:sample_Uniform01(U,S1,S2), (U<P->X=1;X=0).
   59bernoulli(P,0,P1,P2) :- !, P2 is (1-P)*P1.
   60bernoulli(P,1,P1,P2) :- P2 is P*P1.
 binomial(+P:float, +N:natural, -X:natural)// is det
Sample X from a binomial distribution, ie the number of successful trials out of N trials where the probability of success of each trial is P.
   67wrap_rs(3,binomial,sample_Binomial).
   68wrap_disc(3,binomial,prob_Binomial).
 poisson(+A:nonneg, -X:float)// is det
Sample from Poisson distribution of rate A.
   73wrap_rs(2,poisson,sample_Poisson).
   74wrap_disc(2,poisson,prob_Poisson).
 discrete(+A:list(prob), -X:natural)// is det
Sample from a discrete distribution over natural numbers.
   78discrete(Ps,I,p(P1),p(P2))   :- !, nth1(I,Ps,P), P2 is P1*P.
   79discrete(Ps,I,rs(S1),rs(S2)) :- length(Ps,N), plrand:sample_Discrete(N,Ps,I,S1,S2).
 discrete(+O:list(T), +A:list(prob), -X:T)// is det
Sample from a discrete distribution over list of objects.
   84discrete(Xs,Ps,X,rs(S1),rs(S2)) :- !, length(Ps,N), plrand:sample_Discrete(N,Ps,I,S1,S2), nth1(I,Xs,X).
   85discrete(Xs,Ps,X,p(P1),p(P2))   :-
   86   aggregate(sum(P), I^(nth1(I,Xs,X), nth1(I,Ps,P)), Prob),
   87   P2 is P1*Prob.
   88
   89discreteT(Xs,Ps,X,rs(S1),rs(S2)) :- !, functor(Ps,_,N), plrand:sample_DiscreteF(N,Ps,I,S1,S2), arg(I,Xs,X).
   90discreteT(Xs,Ps,X,p(P1),p(P2))   :-
   91   aggregate(sum(P), I^(arg(I,Xs,X), arg(I,Ps,P)), Prob),
   92   P2 is P1*Prob.
 uniform01(-X:float)// is det
Sample X from uniform distribution on [0,1).
   97uniform01(_,p(_),p(0)) :- !.
   98uniform01(_,pd(P),pd(P)) :- !.
   99wrap_rs(1,uniform01,sample_Uniform01).
 normal(-X:float)// is det
Sample from zero-mean unit-variance Gaussian.
  103normal(_,p(_),p(0)) :- !.
  104wrap_rs(1,normal,sample_Normal).
  105wrap_cont(1,normal,prob_Normal).
 exponential(-X:float)// is det
Sample from unit-mean exponential distribution.
  109wrap_rs(1,exponential,sample_Exponential).
  110wrap_cont(1,exponential,prob_Exponential).
 stable(+A, +B, -X:float)// is det
Sample from a Levy-stable distribution.
  115wrap_rs(3,stable,sample_Stable).
 dirichlet(+A:list(nonneg), -X:list(prob))// is det
Sample from a Dirichlet distribution.
  119dirichlet(A,X,rs(S1),rs(S2)) :- !, length(A,N), plrand:sample_Dirichlet(N,A,X,S1,S2).
  120dirichlet(A,X,p(P1),p(P2))   :- length(A,N), plrand:prob_Dirichlet(N,A,X,P), P2 is P1*P.
 uniform(+Items:list(A), -A)// is det
Uniform distribution over a finite number of items. uniform :: list(A) -> expr(A).
  128uniform(O,X,rs(S1),rs(S2)) :- !,
  129   length(O,N),
  130   plrand:sample_Uniform01(U,S1,S2),
  131   I is 1+floor(N*U), nth1(I,O,X).
  132uniform(O,X,p(P1),p(P2)) :-
  133   length(O,N),
  134   aggregate(count,member(X,O),K),
  135   P2 is P1*K/N.
  136
  137uniformT(O,X,p(P1),p(P2)) :- !,
  138   functor(O,_,N),
  139   aggregate(count,I^arg(I,O,X),K),
  140   P2 is K*P1/N.
  141
  142uniformT(O,X,rs(S1),rs(S2)) :-
  143   functor(O,_,N),
  144   plrand:sample_Uniform01(U,S1,S2),
  145   I is 1+floor(N*U), arg(I,O,X).
 uniformP(+P:dcg(-A), -A)// is det
Sample uniformly from all solutions to call(P,X).
  149:- meta_predicate uniformP(3,-,+,-).  150uniformP(P,X) -->
  151   {findall(Y,call(P,Y),YY)},
  152   uniform(YY,X).
 beta(+A:nonneg, +B:nonneg, -X:prob)// is det
Sample from beta distribution.
  156wrap_rs(3,beta,sample_Beta).
  157wrap_cont(3,beta,prob_Beta).
 zeta(+A:nonneg, -X:natural)// is det
Sample from zeta (hyperbolic or power law) distribution over natural numbers. NB: Must have A > 1.
  162wrap_rs(2,zeta,sample_Zeta).
  163wrap_disc(2,zeta,prob_Zeta).
  164% zeta(A,X,)      --> {A>1}, sample_Zeta(A,X1), { X is integer(X1) }.
 gamma(+A:nonneg, -X:float)// is det
Sample from gamma distribution with parameter A.
  168wrap_rs(2,gamma,sample_Gamma).
  169wrap_cont(2,gamma,prob_Gamma).
  170
  171
  172% ^ above use plrand samplers and need randstate
  173% ---------------------- DERIVED DISTRIBUTIONS ---------------------
  174% V below do not use state directly.
 gaussian(+Mean:float, +Var:nonneg, -X:float)// is det
gaussian :: \(float, nonneg) -> expr(float). Sample from Gaussian with given mean and variance.
  179gaussian( Mean, Var, X) --> normal(U), {X is Mean + Var*U}.
 inv_gamma(+A:nonneg, -X:float)// is det
Sample from inverse gamma distribution with parameter A.
  183inv_gamma(A,X)  --> gamma(A,Y), {X is 1/Y}.
 pareto(+A:nonneg, -X:float)// is det
Sample from pareto (power-law) distribution over non-negative reals.
  187pareto(A,X)    --> uniform01(Y), { X is (1-Y)**(-1/A) }.
 students_t(+V:nonneg, -X:float)// is det
Sample from student's t distribution with V degrees of freedom.
  191students_t(V,X)--> {V1 is V/2}, normal(Z), gamma(V1,Y), {X is Z*sqrt(V1/Y)}.
 product_pair(+F:dist(A), +G:dist(B), -X:pair(A,B))// is det
Sample a pair from two independent distributions.
  196product_pair(F,G,X-Y) --> call(F,X), call(G,Y).
 mixture(+Sources:list(expr(A)), +Probs:list(prob), -X:A)// is det
Sample from discrete distribution over Sources with probabilities Probs and then sample from the resulting distribution.

mixture :: \(list(expr(A)), list(prob)) -> expr(A).

  205mixture( Sources, Dist, X) -->
  206   discrete(Dist,I),
  207   {nth1(Sources,I,S)},
  208   call(S,X).
  209
  210prob(Expr,Val,Prob) :- call(Expr,Val,p(1),p(Prob)).
  211pdf(Expr,Val,Prob) :- call(Expr,Val,pd(1),pd(Prob)).
  212prob(Goal,Prob) :- call_dcg(Goal,p(1),p(Prob)).
  213pdf(Goal,Prob)  :- call_dcg(Goal,pd(1),pd(Prob))