/* 

  Hypergeometric 1 distribution in Picat.

  https://en.wikipedia.org/wiki/Hypergeometric_distribution
  """
  (T)he probability of k successes (random draws for which the object 
  drawn has a specified feature) in n draws, without replacement, from 
  a finite population of size N that contains exactly K objects with 
  that feature, wherein each draw is either a success or a failure. 
  In contrast, the binomial distribution describes the probability of 
  k successes in n draws with replacement. 
  """

  Cf https://github.com/distributions-io/hypergeometric-random/blob/master/lib/number.js

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

/*

  Example from
  https://en.wikipedia.org/wiki/Hypergeometric_distribution#Working_example
  """
  Now, assume (for example) that there are 5 green and 45 red marbles in 
  the urn. Standing next to the urn, you close your eyes and draw 10 marbles 
  without replacement. What is the probability that exactly 4 of the 10 
  are green? 
  
  The probability of drawing exactly 4 marbels is: ~0.00396458305801506542
  The probability of drawing exactly 5 marbels is: ~0.00011893749174045196
  """

  This model (method:enumerate) give the following result:
  The probability of drawing exactly 4 marbels is: 30/7567 (0.003964583058015066)
  The probability of drawing exactly 5 marbels is: 9/75670 (0.00011893749174045196

  Which is quite exact...

  Using PDF, CDF, and quantiles:

  Picat> X=hypergeometric1_dist_pdf(50,5,10,4)              
  X = 0.003964583058015

  Picat> X=hypergeometric1_dist_cdf(50,5,10,4)
  X = 0.99988106250826

  Picat> X=hypergeometric1_dist_quantile(50,5,10,0.95)  
  X = 2

  var : c
  Probabilities:
  1: 0.4236000000000000
  0: 0.3169000000000000
  2: 0.2083000000000000
  3: 0.0466000000000000
  4: 0.0044000000000000
  5: 0.0002000000000000
  mean = 0.9986

  var : g
  Probabilities:
  false: 0.9995000000000001
  true: 0.0005000000000000
  mean = [false = 0.9995,true = 0.0005]



*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean]),
  nl,
  % show_store_lengths,
  % fail,
  nl.
go => true.
  
  
model() =>
  K = 5,   % total green marbles: 4 drawn + 1 not drawn
  N = 50,  % total marbles: 5 green + 45 red marbles
  Kk = 4,  % drawn 4 green_marbles
  % Kk 5,  % drawn 5 green_marbles
  Nn = 10, % total drawn green + red marbles
  
  G = hypergeometric1_dist_bool(Kk,N,K,Nn),
  C = hypergeometric1_dist(Kk,N,K,Nn),
  
  add("g",G),
  add("c",C).



/*
   Generate 20 examples from hypergeometric1_dist

   var : cs
   Probabilities (truncated):
   [5,0,1,1,0,0,0,1,0,1,2,3,1,0,0,1,1,1,0,2]: 0.0010000000000000
   [4,2,2,1,1,1,0,4,0,2,1,1,0,2,2,0,0,0,1,0]: 0.0010000000000000
   [4,1,1,1,0,1,0,0,2,0,1,0,0,2,1,3,1,2,1,2]: 0.0010000000000000
   [4,1,1,0,2,1,2,1,1,0,1,1,1,0,0,1,1,1,2,3]: 0.0010000000000000
   .........
   [0,0,0,0,0,1,1,2,0,1,1,1,2,1,0,1,0,0,2,1]: 0.0010000000000000
   [0,0,0,0,0,1,1,1,0,1,1,3,2,1,1,1,1,0,1,1]: 0.0010000000000000
   [0,0,0,0,0,1,0,0,0,2,1,0,0,1,2,2,1,1,2,2]: 0.0010000000000000
   [0,0,0,0,0,0,1,1,1,0,1,1,2,0,2,0,1,1,2,0]: 0.0010000000000000


*/
go2 ?=>
  reset_store,
  run_model(1_000,$model2,[show_probs_trunc]),
  nl,
  % show_store_lengths,
  % fail,
  nl.
go2 => true.
  
  
model2() =>
  K = 5,   % total green marbles: 4 drawn + 1 not drawn
  N = 50,  % total marbles: 5 green + 45 red marbles
  Kk = 4,  % drawn 4 green_marbles
  % Kk 5,  % drawn 5 green_marbles
  Nn = 10, % total drawn green + red marbles
  
  Cs = hypergeometric1_dist_n(Kk,N,K,Nn,20),
  add("cs",Cs).