/* 

  Mixture of gaussian in Picat.

  From BLOG examples/mixture-of-gaussian-full.blog

  Cf my Gamble model gamble_mixture_of_gaussian.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.
% import ordset.

main => go.

/*
  Num accepted samples: 998  Total samples: 144227  (0.007%)
  Num accepted samples: 999  Total samples: 144384  (0.007%)
  Num accepted samples: 1000  Total samples: 144386  (0.007%)
  var : a
  Probabilities (truncated):
  1.844926063142296: 0.0010000000000000
  1.712575484924763: 0.0010000000000000
  1.626403154436083: 0.0010000000000000
  1.549354122350322: 0.0010000000000000
  .........
  -3.528094583343807: 0.0010000000000000
  -3.67080177750047: 0.0010000000000000
  -3.921460785139137: 0.0010000000000000
  -4.000368386535003: 0.0010000000000000
  mean = -0.518013

  var : b
  Probabilities (truncated):
  1.925090707939175: 0.0010000000000000
  1.681389301373006: 0.0010000000000000
  1.594107216372986: 0.0010000000000000
  1.327609607295574: 0.0010000000000000
  .........
  -3.037103253540915: 0.0010000000000000
  -3.400583871047616: 0.0010000000000000
  -3.870876664299051: 0.0010000000000000
  -3.941722411529512: 0.0010000000000000
  mean = 0.0175323

  var : p
  Probabilities (truncated):
  0.995536494810632: 0.0010000000000000
  0.994930161937727: 0.0010000000000000
  0.993423644474295: 0.0010000000000000
  0.993388403459111: 0.0010000000000000
  .........
  0.000010721720421: 0.0010000000000000
  0.000005892033455: 0.0010000000000000
  0.000000357500251: 0.0010000000000000
  0.000000261474443: 0.0010000000000000
  mean = 0.305416

  var : min a,b
  Probabilities (truncated):
  1.228440700943618: 0.0010000000000000
  0.888788127780431: 0.0010000000000000
  0.789545504451513: 0.0010000000000000
  0.743149789452825: 0.0010000000000000
  .........
  -3.870876664299051: 0.0010000000000000
  -3.921460785139137: 0.0010000000000000
  -3.941722411529512: 0.0010000000000000
  -4.000368386535003: 0.0010000000000000
  mean = -0.859597


*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean,
                           % show_percentiles,show_histogram,
                           % show_hpd_intervals,hpd_intervals=[0.84],
                           min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  N = 5,
  P = beta_dist(0.5,1),

  Z = bern_n(P,N),

  A = normal_dist(-1,1),
  B = normal_dist(-1,1),  

  Y = [ cond(Z[I]==1,normal_dist(A,1),normal_dist(B,1)) : I in 1..N],

  Data = [0.2,0.5,1.0,0.5,0.6],
  observe_abc(Data,Y),
  
  if observed_ok then
    add("a",A),
    add("b",B),
    add("p",P),
    add("min a,b",min(A,B)),    
  end.