/* 

  Infinite Gaussian mixtue in Picat.

  From Yi Wu, Lei Li, Stuart Russell, Rastislav Bodik
  "Swift: Compiled Inference for Probabilistic Programming Languages"
  Page 3


  Cf my Gamble model gamble_infinite_gaussian_mixture.rkt

  This program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Picat page: http://www.hakank.org/picat/

*/

import ppl_distributions, ppl_utils.
import util.

main => go.

/*
  var : cluster
  Probabilities:
  3: 0.2410887880751782
  2: 0.2391445236552171
  4: 0.1639662994167207
  1: 0.1471160077770577
  5: 0.1173039533376539
  6: 0.0576798444588464
  7: 0.0207388204795852
  8: 0.0071289695398574
  9: 0.0051847051198963
  11: 0.0006480881399870
  mean = 3.19313

  var : x1
  Probabilities (truncated):
  0.109994485435606: 0.0006480881399870
  0.10998010795709: 0.0006480881399870
  0.109926383586623: 0.0006480881399870
  0.109903149632689: 0.0006480881399870
  .........
  0.090025102284492: 0.0006480881399870
  0.090016177839284: 0.0006480881399870
  0.09001173364283: 0.0006480881399870
  0.090003640334393: 0.0006480881399870
  mean = 0.0994451

  var : x2
  Probabilities (truncated):
  8.197584536865602: 0.0006480881399870
  8.071822350934616: 0.0006480881399870
  7.961080437341861: 0.0006480881399870
  6.918566956916693: 0.0006480881399870
  .........
  -0.254616922842704: 0.0006480881399870
  -0.281961331360441: 0.0006480881399870
  -0.285511979965967: 0.0006480881399870
  -0.300196396925173: 0.0006480881399870
  mean = 1.08385

  var : x3
  Probabilities (truncated):
  8.108589161458903: 0.0006480881399870
  7.959595206214723: 0.0006480881399870
  7.079633964422934: 0.0006480881399870
  6.963562016341839: 0.0006480881399870
  .........
  -0.230320801160971: 0.0006480881399870
  -0.236851737419524: 0.0006480881399870
  -0.273538872759695: 0.0006480881399870
  -0.291008179327642: 0.0006480881399870
  mean = 1.1051


*/
go ?=>
  reset_store,
  run_model(100_000,$model,[show_probs_trunc,mean]),
  nl,
  show_store_lengths,
  % fail,
  nl.
go => true.
  
  
model() =>
  N = 3,
  Cluster = poisson_dist(3)+1, % Number of clusters
  % Mu = normal_dist_n(0,10,Cluster),

  Z = random_integer_n(Cluster,N),

  X = [normal_dist(Z[D],0.1) : D in 1..N],

  observe(abs(X[1]-0.1) <= 0.01), % omit other data points
  % observe(abs(X[2]-0.1) <= 0.01), % omit other data points
  % observe(abs(X[3]-0.1) <= 0.01), % omit other data points  
  
  if observed_ok then
    add("cluster",Cluster),
    add("x1",X[1]),
    add("x2",X[2]), 
    add("x3",X[3]),   
  end.