1/*
    2Seven scientists, posterior estimation in Bayesian models.
    3Suppose seven scientists all go and perform the same experiment, each collecting
    4a measurement xi for i=1,..,7.
    5These scientists are varyingly good at their job, and while we can assume each 
    6scientist would estimate x correctly on average, some of them may have much more
    7error in their measurements than others.
    8They come back with the following seven observations:
    9[-27.020 3.570 8.191 9.898 9.603 9.945 10.056]
   10To model this situation, we put a prior on the mean and the standard deviation
   11of the measurements each of the 7 scientists.
   12For the mean, we use a Gaussian prior with mean 0 and variance 50^2.
   13For the standard deviation, we use a uniform prior between 0 and 25.
   14Given the above measurements, what is the posterior distribution of x?
   15What distribution over noise levels do we infer for each of these scientists' 
   16estimates?
   17From
   18http://www.robots.ox.ac.uk/~fwood/anglican/examples/viewer/?worksheet=gaussian-posteriors
   19*/
   20:- use_module(library(mcintyre)).   21
   22:- if(current_predicate(use_rendering/1)).   23:- use_rendering(c3).   24:- endif.   25:- mc.   26:- begin_lpad.   27
   28val(I,X) :- 
   29  standard_dev(I,Sigma),
   30  mean(M),
   31  measurement(I,M,Sigma,X).
   32% for scientist I we see X sampled from a Gaussian with mean M and standard deviation
   33% Sigma
   34
   35standard_dev(_,S): uniform_dens(S,0,25).
   36% the standard deviation is sampled for all scientists between 0 and 25 
   37% uniformly
   38
   39mean(M): gaussian(M,0, 2500).
   40% the mean is sampled from a Gaussian with mean 0 and variance 50^2
   41
   42measurement(_,M,Sigma,X): gaussian(X,M,Sigma*Sigma).
   43% the measurement is sampled from a Gaussian with mean M and variance
   44% Sigma^2
   45
   46:- end_lpad.   47
   48hist_uncond(Samples,NBins,Chart):-
   49  mc_sample_arg(val(0,X),Samples,X,L0),
   50  histogram(L0,Chart,[nbins(NBins)]).
   51% take Samples samples of X for index 0 (X in val(0,X) and draw a
   52% histogram of the distribution with NBins bins
   53
   54dens_lw(Samples,NBins,Chart,E):-
   55  mc_sample_arg(val(0,Y),Samples,Y,L0),
   56  mc_lw_sample_arg(mean(X),(val(1,-27.020),val(2,3.570),
   57  val(3,8.191),val(4,9.898),val(5,9.603),val(6,9.945),
   58  val(7,10.056)),Samples,X,L),
   59  exp(L,Samples,E),
   60  densities(L0,L,Chart,[nbins(NBins)]).
   61% take Samples samples of X for index 0 (X in val(0,X)) before and after
   62% having observed the scientists' measurements and draw curves of the 
   63% densities using NBins bins
   64
   65chart_lw_noise(Samples,Chart,E):-
   66  mc_lw_sample_arg((standard_dev(1,Y1),standard_dev(2,Y2),standard_dev(3,Y3),standard_dev(4,Y4),
   67    standard_dev(5,Y5),standard_dev(6,Y6),standard_dev(7,Y7)),(val(1,-27.020),val(2,3.570),
   68  val(3,8.191),val(4,9.898),val(5,9.603),val(6,9.945),
   69  val(7,10.056)),Samples,(Y1,Y2,Y3,Y4,Y5,Y6,Y7),L),
   70  exp_noise(L,Samples,E),
   71  E = (E1,E2,E3,E4,E5,E6,E7),
   72  Chart = c3{data:_{x:x, rows:[x-e,1-E1,2-E2,3-E3,4-E4,5-E5,6-E6,7-E7],
   73                    type: bar}}.
   74% take Samples samples of the standard deviation of the measurements of the
   75% scientists (Y1,...,Y7 in standard_dev(1,Y1),...,standard_dev(7,Y7))
   76% given the 7 observations and draw a bar chart of the mean of the samples
   77% for each scientist
   78
   79
   80exp_noise(L,S,(E1,E2,E3,E4,E5,E6,E7)):-
   81  foldl(agg_noise,L,(0,0,0,0,0,0,0),(S1,S2,S3,S4,S5,S6,S7)),
   82  E1 is S1/S,
   83  E2 is S2/S,
   84  E3 is S3/S,
   85  E4 is S4/S,
   86  E5 is S5/S,
   87  E6 is S6/S,
   88  E7 is S7/S.
   89
   90agg_noise((V1,V2,V3,V4,V5,V6,V7)-W,(S1,S2,S3,S4,S5,S6,S7),
   91  (S1+V1*W,S2+V2*W,S3+V3*W,S4+V4*W,S5+V5*W,S6+V6*W,S7+V7*W)).
   92
   93exp(L,S,E):-
   94  foldl(agg,L,0,Sum),
   95  E is Sum/S.
   96
   97agg(V-W,S,S+V*W).