/* 

  Ball selection game in Picat.

  Thanushan Sivapatham:
  "Probability Analysis: Resolving Conditional and Total Probabilities in Tim and David’s Ball Selection Game"
  https://medium.com/puzzle-sphere/probability-analysis-tim-and-davids-ball-selection-game-solving-conditional-and-total-cc6ae1c5cab1
  """
  Friends, Tim and David, are playing a game at a party. 
  Box C contains 3 red balls and 4 black balls, while 
  Box D contains 4 red balls and 3 back balls.
  The ball in boxes A [C!] amd B [D!] are identical in all
  respects, except for their color.
  Tim takes a ball at random from Box C and puts it into
  Box D. Then, David takes a ball at random from Box D.

  Q1) What is the probability that the ball taken out of Box D
      by David is black?
  Q2) What is the probability that the ball taken out of Box C 
      by Tim is black, given that the ball taken out of Box D
      by David is red?
  
  """



  Cf my Gamble model gamble_ball_selection_game.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.

/*
  Q1) What is the probability that the ball taken out of Box D
      by David is black?

      q = 1
      var : draw1
      Probabilities:
      black: 0.5712333333333334
      red: 0.4287666666666667

      var : draw2
      Probabilities:
      red: 0.5548333333333333
      black: 0.4451666666666667


  Q2) What is the probability that the ball taken out of Box C 
      by Tim is black, given that the ball taken out of Box D
      by David is red?

      var : draw1
      Probabilities:
      black: 0.5165987535953979
      red: 0.4834012464046021

      var : draw2
      Probabilities:
      red: 1.0000000000000000


*/
go ?=>
  member(Q,1..2),
  println(q=Q),
  reset_store,
  run_model(30_000,$model(Q),[show_probs_trunc % ,mean % ,
                          %  show_simple_stats % ,
                           % show_percentiles,
                           % show_hpd_intervals,hpd_intervals=[0.84,0.9,0.94,0.99,0.99999],
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  fail,
  nl.
go => true.
  
  
model(Q) =>
  BoxC = [red,red,red,black,black,black,black],
  BoxD = [red,red,red,red,black,black,black],

  Draw1 = uniform_draw(BoxC),
  BoxD2 = BoxD ++ [Draw1],

  Draw2 = uniform_draw(BoxD2),

  if Q == 2 then
    observe(Draw2 == red)
  end,
  
  if Q== 1 ; (Q == 2,observed_ok) then
    add("draw1",Draw1),
    add("draw2",Draw2),    
  end.