/* 

  Hypergeometric 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. 
  """

  hypergeometric_dist mirrors Mathematica's version of Hypergeometric[n, n_succ, n_tot]

  """
  A hypergeometric distribution gives the distribution of the number of successes in 
  n draws from a population of size n_tot containing n_succ successes.
  """

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

/*
  Suppose an urn has 100 elements, of which 40 are special and we draw 10 elements.
  What is the probability of getting 0..10 special elements?
  What is the probability that we get >= 7 special elements?

  Mathematica
  Table((n, PDF(HypergeometricDistribution(10, 40, 100), n)), (n, 1, 10)) 
  -> 
  ((0., 0.00435544),
   (1., 0.0341603), 
   (2., 0.115291), 
   (3., 0.220431), 
   (4., 0.264313), 
   (5., 0.207606), 
   (6., 0.108128), 
   (7., 0.0368556), 
   (8., 0.0078636), 
   (9., 0.000947778), 
   (10., 0.0000489685))

   Probability(x >= 7, x ~ HypergeometricDistribution(10, 40, 100)) ;; N
   -> 0.04571

  
  Picat> X=[K=hypergeometric_dist_pdf(10,40,100,K) : K in 0..10]
  X = [0 = 0.004355440771046,1 = 0.034160319772911,2 = 0.115291079233572,3 = 0.220430742685574,
       4 = 0.264312788683168,5 = 0.207605681292958,6 = 0.108127959006752,7 = 0.036855645175235,
       8 = 0.007863596707647,9 = 0.000947778134255,10 = 0.000048968536937]

var : h
Probabilities:
4: 0.2611000000000000
3: 0.2318000000000000
5: 0.2060000000000000
2: 0.1117000000000000
6: 0.1066000000000000
7: 0.0365000000000000
1: 0.0331000000000000
8: 0.0076000000000000
0: 0.0043000000000000
9: 0.0012000000000000
10: 0.0001000000000000
mean = 3.994
Percentiles:
(0.001 0)
(0.01 1)
(0.025 1)
(0.05 2)
(0.1 2)
(0.25 3)
(0.5 4)
(0.75 5)
(0.84 5)
(0.9 6)
(0.95 6)
(0.975 7)
(0.99 7)
(0.999 9)
(0.9999 9.0000999999992928)
(0.99999 9.9000100000012026)
HPD intervals:
HPD interval (0.5): 2.00000000000000..4.00000000000000
HPD interval (0.84): 1.00000000000000..5.00000000000000
HPD interval (0.9): 2.00000000000000..6.00000000000000
HPD interval (0.95): 1.00000000000000..6.00000000000000
HPD interval (0.99): 0.00000000000000..7.00000000000000
HPD interval (0.999): 0.00000000000000..9.00000000000000
HPD interval (0.9999): 0.00000000000000..9.00000000000000
HPD interval (0.99999): 0.00000000000000..10.00000000000000
Histogram:
      0:    43 #                                                            (0.004 / 0.004)
      1:   331 ########                                                     (0.033 / 0.037)
      2:  1117 ##########################                                   (0.112 / 0.149)
      3:  2318 #####################################################        (0.232 / 0.381)
      4:  2611 ############################################################ (0.261 / 0.642)
      5:  2060 ###############################################              (0.206 / 0.848)
      6:  1066 ########################                                     (0.107 / 0.955)
      7:   365 ########                                                     (0.036 / 0.991)
      8:    76 ##                                                           (0.008 / 0.999)
      9:    12                                                              (0.001 / 1.000)
     10:     1                                                              (0.000 / 1.000)

var : p
Probabilities:
false: 0.9546000000000000
true: 0.0454000000000000
mean = [false = 0.9546,true = 0.0454]
Percentiles:
(0.001 false)
(0.01 false)
(0.025 false)
(0.05 false)
(0.1 false)
(0.25 false)
(0.5 false)
(0.75 false)
(0.84 false)
(0.9 false)
(0.95 false)
(0.975 true)
(0.99 true)
(0.999 true)
(0.9999 true)
(0.99999 true)
HPD intervals:
show_hpd_intervals: data is not numeric
Histogram:
false     :  9546 ############################################################ (0.955 / 0.955)
true      :   454 ###                                                          (0.045 / 1.000)


*/
go ?=>
  reset_store,
  run_model(10_000,$model),
  nl,
  % show_store_lengths,
  % fail,
  nl.
go => true.
  
model() =>
  H = hypergeometric_dist(10,40,100),
  P = check(H >= 7),
  add("h",H),
  add("p",P).