```    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
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)```