/* 

  Two aces in Picat.

  From Statistics101 (Resample Stats)
  http://www.statistics101.net/QuickReference.pdf
  Page 20
  """
  What is the probability of drawing an ace at random from a
  deck of cards, and then on your second draw drawing another ace?
  ...
  Thus the answer to the question would be 4/52 x 3/51 = .00452.
  """


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

/*
  Using resampling:

  var : p
  Probabilities:
  false: 0.9958000000000000
  true: 0.0042000000000000
  mean = 0.0042

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean % ,
                           % show_percentiles,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  Deck = rep(4,1..13).flatten,
  Sample = draw_without_replacement(2,Deck),

  P = check(Sample == [1,1]),
  
  add("p",P).


/*
  Using hypergeometric distribution

  var : p1
  Probabilities:
  0.00452489: 1.0000000000000000
  mean = 0.00452489

  var : p2
  Probabilities:
  false: 0.9953000000000000
  true: 0.0047000000000000
  mean = 0.0047


  Picat> P1 = hypergeometric_dist_pdf(2,4,52,2) 
  P1 = 0.004524886877828

  Picat> pdf_all($hypergeometric_dist(2,4,52),0.01,0.999).printf_list 
   0 0.850678733031662 
   1 0.144796380090491 
   2 0.004524886877828 

*/
go2 ?=>
  reset_store,
  run_model(10_000,$model2,[show_probs_trunc,mean % ,
                           % show_percentiles,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go2 => true.
  
  
model2() =>
  P1 = hypergeometric_dist_pdf(2,4,52,2),
  P2 = check(2 == hypergeometric_dist(2,4,52)),  
  
  add("p1",P1),
  add("p2",P2).