1/*
    2Mixture model from a Dirichlet process (DP), see http://www.robots.ox.ac.uk/~fwood/anglican/examples/viewer/?worksheet=nonparametrics/dp-mixture-model
    3https://en.wikipedia.org/wiki/Dirichlet_process
    4Samples are drawn from a mixture of normal distributions whose parameters are
    5defined by means of a Dirichlet process, so the number of components is not
    6fixed in advance. For each component, the variance is sampled from a gamma
    7disrtibution and the mean is sampled from a Guassian with mean 0 and variance
    830 times the variance of the compoment.
    9Given some observations, the aim is to find how the distribution of values is
   10updated. Less observations are considered with respect to http://www.robots.ox.ac.uk/~fwood/anglican/examples/viewer/?worksheet=nonparametrics/dp-mixture-model
   11because the weights go rapidly to 0.
   12*/

?- prior(200,100,G). % draw the prior density ?- post(200,100,G). % draw the posterior density

*/

   20 :- use_module(library(mcintyre)).   21:- if(current_predicate(use_rendering/1)).   22:- use_rendering(c3).   23:- endif.   24:- mc.   25:- begin_lpad.   26
   27dp_n_values(N,N,_Alpha,[]):-!.
   28
   29dp_n_values(N0,N,Alpha,[[V]-1|Vs]):-
   30  N0<N,
   31  dp_value(N0,Alpha,V),
   32  N1 is N0+1,
   33  dp_n_values(N1,N,Alpha,Vs).
   34  
   35dp_value(NV,Alpha,V):-
   36  dp_stick_index(NV,Alpha,I),
   37  dp_pick_value(I,NV,V).
   38
   39dp_pick_value(I,NV,V):-
   40  ivar(I,IV),
   41  Var is 1.0/IV,
   42  mean(I,Var,M),
   43  value(NV,M,Var,V).
   44
   45ivar(_,IV):gamma(IV,1,0.1).
   46
   47mean(_,V0,M):gaussian(M,0,V):-V is V0*30.
   48
   49value(_,M,V,Val):gaussian(Val,M,V).
   50
   51dp_stick_index(NV,Alpha,I):-
   52  dp_stick_index(1,NV,Alpha,I).
   53
   54dp_stick_index(N,NV,Alpha,V):-
   55  stick_proportion(N,Alpha,P),
   56  choose_prop(N,NV,Alpha,P,V).
   57  
   58choose_prop(N,NV,_Alpha,P,N):-
   59  pick_portion(N,NV,P).
   60
   61choose_prop(N,NV,Alpha,P,V):-
   62  neg_pick_portion(N,NV,P),
   63  N1 is N+1,
   64  dp_stick_index(N1,NV,Alpha,V).
   65 
   66
   67
   68stick_proportion(_,Alpha,P):beta(P,1,Alpha).
   69
   70pick_portion(N,NV,P):P;neg_pick_portion(N,NV,P):1-P.
   71
   72:- end_lpad.   73
   74
   75obs([-1,7,3]).
   76
   77prior(Samples,NBins,Chart):-
   78  mc_sample_arg_first(dp_n_values(0,Samples,10.0,V),1,V,L),
   79  L=[Vs-_],
   80  histogram(Vs,Chart,[nbins(NBins)]).
   81
   82post(Samples,NBins,Chart):-
   83  obs(O),
   84  maplist(to_val,O,O1),
   85  length(O1,N),
   86  mc_lw_sample_arg_log(dp_value(0,10.0,T),dp_n_values(0,N,10.0,O1),Samples,T,L),
   87  maplist(keys,L,LW),
   88  min_list(LW,Min),
   89  maplist(exp(Min),L,L1),
   90  density(L1,Chart,[nbins(NBins),min(-8),max(15)]).
   91
   92keys(_-W,W).
   93
   94exp(Min,L-W,L-W1):- W1 is exp(W-Min).
   95
   96to_val(V,[V]-1)