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
41   length(A1, Arity),
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),
48   append(A1, [P],           PredArgs), Body =.. [Pred | PredArgs].
49
51   length(A1, Arity),
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))