1/* 2Posterior estimation in Bayesian models. 3We are trying to estimate the true value of a Gaussian distributed random 4variable, given some observed data. The variance is known (2) and we 5suppose that the mean has a Gaussian distribution with mean 1 and variance 65. We take different measurement (e.g. at different times), indexed 7with an integer. 8Given that we observe 9 and 8 at indexes 1 and 2, how does the distribution 9of the random variable (value at index 0) changes with respect to the case of 10no observations? 11From 12http://www.robots.ox.ac.uk/~fwood/anglican/examples/viewer/?worksheet=gaussian-posteriors 13*/ 14:- use_module(library(mcintyre)). 15 16:- if(current_predicate(use_rendering/1)). 17:- use_rendering(c3). 18:- endif. 19 20 21 22 23:- mc. 24:- begin_lpad. 25 26val(I,X) :- 27 mean(M), 28 val(I,M,X). 29% at time I we see X sampled from a Gaussian with mean M and variamce 2.0 30 31mean(M) user(M,gauss(1.0, 5.0)). 32% Gaussian distribution of the mean of the Gaussian of the variable 33 34val(_,M,X) user(X,gauss(M, 2.0)). 35% Gaussian distribution of the variable 36 37 38:- end_lpad. 39 40gauss(Mean,Variance,S):- 41 number(Mean),!, 42 random(U1), 43 random(U2), 44 R is sqrt(-2*log(U1)), 45 Theta is 2*pi*U2, 46 S0 is R*cos(Theta), 47 StdDev is sqrt(Variance), 48 S is StdDev*S0+Mean. 49 50gauss(Mean,Variance,S,D):- 51 StdDev is sqrt(Variance), 52 D is 1/(StdDev*sqrt(2*pi))*exp(-(S-Mean)*(S-Mean)/(2*Variance)). 53 54hist_uncond(Samples,NBins,Chart):- 55 mc_sample_arg(val(0,X),Samples,X,L0), 56 histogram(L0,Chart,[nbins(NBins)]). 57% plot an histogram of the density of the random variable before any 58% observations by taking Samples samples and by dividing the domain 59% in NBins bins 60 61dens_lw(Samples,NBins,Chart):- 62 mc_sample_arg(val(0,Y),Samples,Y,L0), 63 mc_lw_sample_arg(val(0,X),(val(1,9),val(2,8)),Samples,X,L), 64 densities(L0,L,Chart,[nbins(NBins)]). 65% plot the densities of the random variable before and after 66% observing 9 and 8 by taking Samples samples using likelihood weighting 67% and by dividing the domain 68% in NBins bins 69 70dens_part(Samples,NBins,Chart):- 71 mc_sample_arg(val(0,Y),Samples,Y,L0), 72 mc_particle_sample_arg(val(0,X),[val(1,9),val(2,8)],Samples,X,L), 73 densities(L0,L,Chart,[nbins(NBins)]). 74% plot the densities of the random variable before and after 75% observing 9 and 8 by taking Samples samples using particle filtering 76% and by dividing the domain 77% in NBins bins
?-
dens_lw(1000,40,G)
. % plot the densities of the random variable before and after % observing 9 and 8 using likelihood weighting ?-dens_part(1000,40,G)
. % plot the densities of the random variable before and after % observing 9 and 8 using particle filtering ?-hist_uncond(10000,40,G)
. % plot an histogram of the density of the random variable before any % observations ?-mc_lw_expectation(val(0,X),(val(1,9),val(2,8)),1000,X,E)
. % E = 7.166960047178755 ?-mc_expectation(val(0,X),10000,X,E)
. % E = 0.9698875384639362.*/