/* 

  Coins learning in Picat.

  From srl-pp-tutorial-wasp-stockholm.pdf
  "Statistical Relational learning and Probabilistic Programming"
  by Luc De Raedt, Anton Dries, and Angelika Kimmig
  https://dtai.cs.kuleuven.be/problog/wasp17-tutorial.html
  
  slide 394f:
  """
  - Flipping a coin with unknown weight
  - Prior: uniform distribution on (0,1)
  - Observation: 5x heads in a row
  """

  The ProbLog model return the following which corresponds to the
  density of the probabilist 
  """
   weight(c1,0.1): 3.3994697e-13
   weight(c1,0.3): 2.1679411e-06
   weight(c1,0.5): 0.0041497433
   weight(c1,0.7): 0.1317485 
   weight(c1,0.9): 0.86409959         <----

   weight(c2,0.1): 3.2276726e-06
   weight(c2,0.3): 0.024109997
   weight(c2,0.5): 0.66724754         <----
   weight(c2,0.7): 0.30628626
   weight(c2,0.9): 0.0023529733
  """

  Here we observe 13 tosses and detects the probability of throwing a head (1).
  The corresponding values of the ProbLog model is in "params ix".

  Note that this model is much simpler than the ProbLog model since ProbLog 
  don't support continous distributions. 
 
  For a model which is closer to the ProbLog model, see ppl_coins_learning2.pi

  Cf 
  - ppl_coins_learning2.pi
  - my Gamble model gamble_coins_learning.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, ppl_common_utils.
import util.
% import ordset.

main => go.

/*

  data = [1,1,1,1,1,1,1,1,1,1,1,1,1]
  var : p
  Probabilities (truncated):
  0.999979366082688: 0.0010000000000000
  0.999904904048846: 0.0010000000000000
  0.999808776657939: 0.0010000000000000
  0.999781499151039: 0.0010000000000000
  .........
  0.653788326612575: 0.0010000000000000
  0.639880772978012: 0.0010000000000000
  0.625929738686387: 0.0010000000000000
  0.601081719436255: 0.0010000000000000
  mean = 0.931989


  data = [1,0,1,1,1,1,1,0,0,1,0,0,1]
  var : p
  Probabilities (truncated):
  0.877377217578412: 0.0010000000000000
  0.860880428860374: 0.0010000000000000
  0.860031684329748: 0.0010000000000000
  0.854907264399765: 0.0010000000000000
  .........
  0.26219767856514: 0.0010000000000000
  0.254154563068484: 0.0010000000000000
  0.224604931764586: 0.0010000000000000
  0.219217389924087: 0.0010000000000000
  mean = 0.601261


*/
go ?=>
  Data1 = [1,1,1,1,1,1,1,1,1,1,1,1,1],
  Data2 = [1,0,1,1,1,1,1,0,0,1,0,0,1],
  member(Data,[Data1,Data2]),
  println(data=Data),
  reset_store,
  run_model(10_000,$model(Data),[show_probs_trunc,mean,
                           min_accepted_samples=1000 % ,show_accepted_samples=true
                          ]),
  nl,
  fail,
  nl.
go => true.
  
  
model(Data) =>
  P = uniform(0,1),

  Y = bern_n(P,Data.len),
  
  observe_abc(Data,Y,1/10),
  
  if observed_ok then
    add("p",P)
  end.