1/* 2 * Prolog part of random generator library 3 * Samer Abdallah (2009) 4*/ 5 6:- module(randmeta, [ 7 rv/2 % -Spec, -Type 8 , sample/4 % +Dist, -Value, +StateIn, -StateOut 9 , sample/2 % +Dist, -Value 10 , dps_dist/3 11 12 , op(900,fx,\\) 13 ]).
21:- use_module(library(plrand)). 22:- multifile sample/4, rv/2.
34rv( raw, natural). 35rv( uniform01, nonneg). 36rv( normal, real). 37rv( exponential, nonneg). 38rv( gamma(nonneg), nonneg). 39rv( students_t(nonneg), real). 40rv( poisson(nonneg), natural). 41rv( invgamma(nonneg), nonneg). 42rv( beta(nonneg,nonneg), nonneg). 43rv( zeta(nonneg), nonneg). 44rv( binomial(nonneg,natural), natural). 45rv( dirichlet(natural,list(nonneg)), list(nonneg)). 46rv( dirichlet(list(nonneg)), list(nonneg)). 47rv( dirichlet(tuple(nonneg)), tuple(nonneg)). 48rv( stable(nonneg,real), real). 49rv( bernoulli(nonneg), atom). 50rv( bernoulli(nonneg), atom). 51rv( discrete(natural,list(nonneg)),natural). 52rv( discrete(list(nonneg)),natural). 53rv( discrete(tuple(nonneg)),natural).
sample/2 uses and modifies the global state. sample/4 uses a given random generator state and returns the state on completion, and is designed to compatible with DCG syntax to hide the threading of the random state through several consecutive calls.
DistExpression is an expression describing a distribution. The head functor can be a distribution name (as listed by rv/2) or one of a number of arithmetic operators or term constructors. The arguments (in almost all cases) can then be further DistExpressions which are evaluated recursively. Valid non-distributional termsare:
return samples from X, Y etc in a list
rep(N,X)returns N independent sample from X in a list (N must be a constant)
factorial(N,X)
returns N independent sample from X in an N argument term \(X1,...,XN)For example
?- sample(invgamma(1)*rep(8,discrete(dirichlet(rep(4,0.5)))),X). X = [3, 2, 3, 3, 3, 1, 2, 2] .
99sample( raw, X) --> plrand:sample_Raw(X1), {X is integer(X1)}. 100 101sample( uniform01, X) --> plrand:sample_Uniform01(X). 102sample( normal, X) --> plrand:sample_Normal(X). 103sample( exponential, X) --> plrand:sample_Exponential(X). 104sample( gamma(A), X) --> sample(A,A1), plrand:sample_Gamma(A1,X). 105sample( poisson(A), X) --> sample(A,A1), plrand:sample_Poisson(A1,X). 106sample( invgamma(A), X) --> sample(A,A1), plrand:sample_Gamma(A1,Y), {X is 1/Y}. 107sample( beta(A,B), X) --> sample(A,A1), sample(B,B1), plrand:sample_Beta(A1,B1,X). 108sample( zeta(A), X) --> sample(A,A1), plrand:sample_Zeta(A1,X1), {X is integer(X1)}. 109sample( pareto(A), X) --> sample(A,A1), plrand:sample_Uniform01(Y), {X is (1-Y)**(-1/A1) }. 110sample( binomial(P,N), X) --> sample(P,P1), sample(N,N1), plrand:sample_Binomial(P1,N1,X). 111 112sample( stable(A,B), X) --> sample(A,A1), sample(B,B1), plrand:sample_Stable(A1,B1,X). 113 114sample( dirichlet(N,A), X) --> sample(A,A1), plrand:sample_Dirichlet(N,A1,X). 115sample( dirichlet(A), X) --> sample(A,A1), 116 ( {A1=[_|_]} 117 -> {length(A1,N)}, plrand:sample_Dirichlet(N,A1,X) 118 ; {functor(A1,F,N), functor(X,F,N)}, plrand:sample_DirichletF(N,A1,X) 119 ). 120 121sample( discrete(N,P), X) --> sample(P,P1), plrand:sample_Discrete(N,P1,X). 122sample( discrete(P), X) --> 123 sample(P,P1), 124 ( {P1=[_|_]} 125 -> {length(P1,N)}, plrand:sample_Discrete(N,P1,X) 126 ; {functor(P1,_,N)}, plrand:sample_DiscreteF(N,P1,X) 127 ). 128 129sample( bernoulli(P), X) --> sample(P,P1), plrand:sample_Uniform01(U), {U<P1->X=1;X=0}. 130sample( students_t(V), X) --> sample(V/2,V1), sample(normal*sqrt(V1/gamma(V1)),X). 131 132% dps(Vals) represents a countably infinite discrete distribution. It is an infinite 133% stream of weight:value pairs. 134 135sample(dps([B1:X1|Vals]),X) --> 136 plrand:sample_Uniform01(U), 137 ( {U<B1} -> {X=X1} 138 ; sample(dps(Vals),X) 139 ). 140 141 142sample( X*Y, Z) --> sample(X,X1), sample(Y,Y1), {Z is X1*Y1}. 143sample( X/Y, Z) --> sample(X,X1), sample(Y,Y1), {Z is X1/Y1}. 144sample( X+Y, Z) --> sample(X,X1), sample(Y,Y1), {Z is X1+Y1}. 145sample( X-Y, Z) --> sample(X,X1), sample(Y,Y1), {Z is X1-Y1}. 146sample( -X, Z) --> sample(X,X1), {Z is -X1}. 147sample( sqrt(X), Z) --> sample(X,X1), {Z is sqrt(X1)}. 148 149sample( [], []) --> []. 150sample( [X|XX], [Z|ZZ]) --> sample(X,Z), sample(XX,ZZ). 151sample( rep(N,X), Z) --> {length(Z,N)}, seqmap(sample(X),Z). 152sample( stream(X), Z) --> spawn(S), {freeze(Z,unfold_stream(S,X,Z))}. 153sample( factorial(N,X), Z) --> {functor(Z,\,N)}, seqmapargs(sample(X),N,Z). 154sample( \\(G), X) --> call(G,X). 155sample( Tuple, Value) --> 156 {functor(Tuple,F,N), functor(Value,F,N), tuple_functor(F)}, 157 seqmapargs(sample,N,Tuple,Value). 158 159sample( N, N, S, S) :- number(N), !. 160 161sample(M,X) :- get_rnd_state(S1), sample(M,X,S1,S2), set_rnd_state(S2). 162 163unfold_stream(S1,X,[Z1|ZX]) :- sample(X,Z1,S1,S2), freeze(ZX,unfold_stream(S2,X,ZX)). 164 165% Return truncated infinite discrete distribution 166dps_dist(dps(L),Probs,Vals) :- unfreeze_dps(1,L,Probs,Vals). 167 168unfreeze_dps(_,_,[],[]). 169unfreeze_dps(P0,[Q:V|T],[P|PX],[V|VX]) :- 170 P is P0*Q, P1 is P0*(1-Q), 171 unfreeze_dps(P1,T,PX,VX). 172 173 174tuple_functor(\). 175tuple_functor(tuple). 176tuple_functor(vec). 177 178seqmapargs(P,N,X1) --> 179 ( {N>0} 180 -> {succ(M,N), arg(N,X1,X1N)}, 181 call(P,X1N), 182 seqmapargs(P,M,X1) 183 ; [] 184 ). 185seqmapargs(P,N,X1,X2) --> 186 ( {N>0} 187 -> {succ(M,N), arg(N,X1,X1N), arg(N,X2,X2N)}, 188 call(P,X1N,X2N), 189 seqmapargs(P,M,X1,X2) 190 ; [] 191 ). 192 193seqmap(_,[]) --> []. 194seqmap(P,[A|AX]) --> call(P,A), seqmap(P,AX). 195seqmap(_,[],[]) --> []. 196seqmap(P,[A|AX],[B|BX]) --> call(P,A,B), seqmap(P,AX,BX)
Random variable sampling interpreter
This module provides a mini language for describing and sampling from random variables. */