1/*
    2Dirichlet process (DP), see https://en.wikipedia.org/wiki/Dirichlet_process
    3Samples are drawn from a base distribution. New samples have a nonzero
    4probability of being equal to already sampled values. The process depends
    5on a parameter alpha (concentration parameter): with alpha->0, a single 
    6value is sampled, with alpha->infinite the distribution is equal to the base
    7distribution.
    8In this example the base distribution is a Gaussian with mean 0 and variance
    91, as in https://en.wikipedia.org/wiki/Dirichlet_process#/media/File:Dirichlet_process_draws.svg
   10To model the process, this example uses a stick breaking process: to sample
   11a value, a sample beta_1 is taken from Beta(1,alpha) and a coin with heads
   12probability beta_1 is flipped. If the coin lands heads, a sample from the base
   13distribution is taken and returned. Otherwise, a sample beta_2 is taken again
   14from Beta(1,alpha) and a coin is flipped. This procedure is repeated until
   15a heads is obtained, the index of i beta_i is the index of the value to be 
   16returned.
   17The example queries show both the distribution of indexes and values of the DP.
   18Moreover, they show the distribution of unique indexes as in 
   19http://www.robots.ox.ac.uk/~fwood/anglican/examples/viewer/?worksheet=nonparametrics/dp-mixture-model
   20*/

   32 :- use_module(library(mcintyre)).   33
   34:- if(current_predicate(use_rendering/1)).   35:- use_rendering(c3).   36:- endif.   37:- mc.   38:- begin_lpad.   39
   40% dp_n_values(N0,N,Alpha,L)
   41% returns in L a list of N-N0 samples from the DP with concentration parameter
   42% Alpha
   43dp_n_values(N,N,_Alpha,[]):-!.
   44
   45dp_n_values(N0,N,Alpha,[[V]-1|Vs]):-
   46  N0<N,
   47  dp_value(N0,Alpha,V),
   48  N1 is N0+1,
   49  dp_n_values(N1,N,Alpha,Vs).
   50  
   51% dp_value(NV,Alpha,V)
   52% returns in V the NVth sample from the DP with concentration parameter
   53% Alpha
   54dp_value(NV,Alpha,V):-
   55  dp_stick_index(NV,Alpha,I),
   56  dp_pick_value(I,V).
   57
   58% dp_pick_value(I,V)
   59% returns in V the value of index I of the base distribution 
   60% (in this case N(0,1))
   61dp_pick_value(_,V):gaussian(V,0,1).
   62
   63% dp_stick_index(NV,Alpha,I)
   64% returns in I the index of the NVth sample from the DP
   65dp_stick_index(NV,Alpha,I):-
   66  dp_stick_index(1,NV,Alpha,I).
   67
   68dp_stick_index(N,NV,Alpha,V):-
   69  stick_proportion(N,Alpha,P),
   70  choose_prop(N,NV,Alpha,P,V).
   71  
   72% choose_prop(N,NV,Alpha,P,V)
   73% returns in V the index of the end of the stick breaking process starting
   74% from index N for the NVth value to be sampled from the DP
   75choose_prop(N,NV,_Alpha,P,N):-
   76  pick_portion(N,NV,P).
   77
   78choose_prop(N,NV,Alpha,P,V):-
   79  neg_pick_portion(N,NV,P),
   80  N1 is N+1,
   81  dp_stick_index(N1,NV,Alpha,V).
   82 
   83% sample of the beta_i parameters
   84stick_proportion(_,Alpha,P):beta(P,1,Alpha).
   85
   86% flip of the coin for the portion of the stick of size P
   87pick_portion(N,NV,P):P;neg_pick_portion(N,NV,P):1-P.
   88
   89:- end_lpad.   90
   91hist(Samples,NBins,Chart):-
   92  mc_sample_arg(dp_stick_index(1,10.0,V),Samples,V,L),
   93  histogram(L,Chart,[nbins(NBins)]).
   94
   95hist_repeated_indexes(Samples,NBins,Chart):-
   96  repeat_sample(0,Samples,L),
   97  histogram(L,Chart,[nbins(NBins)]).
   98
   99repeat_sample(S,S,[]):-!.
  100
  101repeat_sample(S0,S,[[N]-1|LS]):-
  102  mc_sample_arg_first(dp_stick_index(1,1,10.0,V),10,V,L),
  103  length(L,N),
  104  S1 is S0+1,
  105  repeat_sample(S1,S,LS).
  106
  107hist_val(Samples,NBins,Chart):-
  108  mc_sample_arg_first(dp_n_values(0,Samples,10.0,V),1,V,L),
  109  L=[Vs-_],
  110  histogram(Vs,Chart,[nbins(NBins)])