/* 

  Conditional probability in Picat.

  From 
  "<lambda>PSI - Exact Inference for Higher Order Probabilistic Programs"
  https://www.youtube.com/watch?v=GraKkfALsFY
  @0:12

  PSI model
  """
  def main() {
    p := uniform(0,1);
    r := 0;
    for i in [0..10] {
      r += flip(p);
    }

    observe(r % 3 == 0);
  
    return (p);
  
  }
  """

  I.e. what is the probability that a random generated series
  of 10 0/1 is 0 modulo 3?

  Cf my Gamble model gamble_youtube1.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 : p
  Probabilities (truncated):
  0.999516086186988: 0.0000916758342501
  0.998140586539237: 0.0000916758342501
  0.997773592359281: 0.0000916758342501
  0.997011575380811: 0.0000916758342501
  .........
  0.000231524463851: 0.0000916758342501
  0.000202067196463: 0.0000916758342501
  0.000125541351794: 0.0000916758342501
  0.000088965985965: 0.0000916758342501
  mean = 0.455435

  var : r
  Probabilities:
  3: 0.2553171983865053
  6: 0.2519251925192519
  0: 0.2500000000000000
  9: 0.2427576090942428
  mean = 4.46232

*/
go ?=>
  reset_store,
  run_model(30_000,$model,[show_probs_trunc,mean]),
  nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  P = uniform_dist(0,1),
  R = [bernoulli_dist(P) : _ in 1..10].sum,
  
  observe(R mod 3 == 0),
  
  if observed_ok then
    add("p",P),
    add("r",R)
  end.