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 ]).
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].
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.
67wrap_rs(3,binomial,sample_Binomial). 68wrap_disc(3,binomial,prob_Binomial).
73wrap_rs(2,poisson,sample_Poisson). 74wrap_disc(2,poisson,prob_Poisson).
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).
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.
97uniform01(_,p(_),p(0)) :- !. 98uniform01(_,pd(P),pd(P)) :- !. 99wrap_rs(1,uniform01,sample_Uniform01).
103normal(_,p(_),p(0)) :- !. 104wrap_rs(1,normal,sample_Normal). 105wrap_cont(1,normal,prob_Normal).
109wrap_rs(1,exponential,sample_Exponential). 110wrap_cont(1,exponential,prob_Exponential).
115wrap_rs(3,stable,sample_Stable).
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.
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).
call(P,X)
.149:- meta_predicate uniformP( , , , ). 150uniformP(P,X) --> 151 {findall(Y,call(P,Y),YY)}, 152 uniform(YY,X).
156wrap_rs(3,beta,sample_Beta). 157wrap_cont(3,beta,prob_Beta).
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) }.
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.
expr(float)
.
Sample from Gaussian with given mean and variance.
179gaussian( Mean, Var, X) --> normal(U), {X is Mean + Var*U}.
183inv_gamma(A,X) --> gamma(A,Y), {X is 1/Y}.
187pareto(A,X) --> uniform01(Y), { X is (1-Y)**(-1/A) }.
191students_t(V,X)--> {V1 is V/2}, normal(Z), gamma(V1,Y), {X is Z*sqrt(V1/Y)}.
196product_pair(F,G,X-Y) --> call(F,X), call(G,Y).
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))
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.