/* 

  Negative binomial test in Picat.

  From BLOG documentation:
  """
  NegativeBinomial distribution generates the number of failures until the rth success 
  in a sequence of independent Bernoulli trials each with probability of success p.

  Example: The following code defines a random function symbol x distributed according to a
  Negative Binomial distribution with probability of success p = 0.8 and number of failures r = 2.
  random Integer x ~ NegativeBinomial(2, 0.8);
  """


  https://stattrek.com/probability-distributions/negative-binomial.aspx
  """
  Bob is a high school basketball player. He is a 70% free throw shooter. That means 
  his probability of making a free throw is 0.70. During the season, what is the 
  probability that Bob makes his third free throw on his fifth shot?

  Solution: This is an example of a negative binomial experiment. The probability of 
  success (P) is 0.70, the number of trials (x) is 5, and the number of 
  successes (r) is 3.

  To solve this problem, we enter these values into the negative binomial formula.
  
  b*(x; r, P) = x-1Cr-1 * Pr * Qx - r
  b*(5; 3, 0.7) = 4C2 * 0.73 * 0.32
  b*(5; 3, 0.7) = 6 * 0.343 * 0.09 = 0.18522
  
  Thus, the probability that Bob will make his third successful free throw on his 
  fifth shot is 0.18522.
  """

  Cf my Gamble model gamble_negative_binomial_test.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:
false: 0.8112000000000000
true: 0.1888000000000000
mean = [false = 0.8112,true = 0.1888]
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 true)
(0.9 true)
(0.95 true)
(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     :  8112 ############################################################ (0.811 / 0.811)
true      :  1888 ##############                                               (0.189 / 1.000)

var : y
Probabilities:
0: 0.3354000000000000
1: 0.3127000000000000
2: 0.1888000000000000
3: 0.0925000000000000
4: 0.0430000000000000
5: 0.0164000000000000
6: 0.0075000000000000
7: 0.0021000000000000
8: 0.0013000000000000
9: 0.0002000000000000
11: 0.0001000000000000
mean = 1.2948
Percentiles:
(0.001 0)
(0.01 0)
(0.025 0)
(0.05 0)
(0.1 0)
(0.25 0)
(0.5 1)
(0.75 2)
(0.84 3)
(0.9 3)
(0.95 4)
(0.975 5)
(0.99 6)
(0.999 8)
(0.9999 9.0001999999985856)
(0.99999 10.800020000002405)
HPD intervals:
HPD interval (0.5): 0.00000000000000..1.00000000000000
HPD interval (0.84): 0.00000000000000..3.00000000000000
HPD interval (0.9): 0.00000000000000..3.00000000000000
HPD interval (0.95): 0.00000000000000..4.00000000000000
HPD interval (0.99): 0.00000000000000..6.00000000000000
HPD interval (0.999): 0.00000000000000..8.00000000000000
HPD interval (0.9999): 0.00000000000000..9.00000000000000
HPD interval (0.99999): 0.00000000000000..11.00000000000000
Histogram:
      0:  3354 ############################################################ (0.335 / 0.335)
      1:  3127 ########################################################     (0.313 / 0.648)
      2:  1888 ##################################                           (0.189 / 0.837)
      3:   925 #################                                            (0.092 / 0.929)
      4:   430 ########                                                     (0.043 / 0.972)
      5:   164 ###                                                          (0.016 / 0.989)
      6:    75 #                                                            (0.007 / 0.996)
      7:    21                                                              (0.002 / 0.998)
      8:    13                                                              (0.001 / 1.000)
      9:     2                                                              (0.000 / 1.000)
     11:     1                                                              (0.000 / 1.000)

*/
go ?=>
  reset_store,
  run_model(10_000,$model),
  nl,
  % show_store_lengths,
  % fail,
  nl.
go => true.
  
  
model() =>
  Y = negative_binomial_dist(3,0.7),

  P = check(Y == 2),
  
  add("y",Y),
  add("p",P).


/*
  Using 0/1 (bernoulli) representation:

  var : p
  Probabilities:
  false: 0.81579999999999997 (4079 / 5000)
  true: 0.18420000000000000 (921 / 5000)
  mean = [false = 0.8158,true = 0.1842]

*/
go2 ?=>
  reset_store,
  run_model(30_000,$model2,[show_probs_rat,mean]),
  nl,
  % show_store_lengths,
  % fail,
  nl.
go2 => true.
  
  
model2() =>
  N = 5,
  Hits = bern_n(0.7,N),
  P = check((Hits.sum == 3, Hits[N] == 1)),
  add("p",P).