/* 

  Beta binomial distribution in Picat.

  https://en.wikipedia.org/wiki/Beta-binomial_distribution#Generating_beta_binomial-distributed_random_variables
  """
  The beta-binomial distribution is the binomial distribution in which the 
  probability of success at each of n trials is not fixed but randomly drawn 
  from a beta distribution. It is frequently used in Bayesian statistics, 
  empirical Bayes methods and classical statistics to capture overdispersion 
  in binomial type distributed data.

  The beta-binomial is a one-dimensional version of the Dirichlet-multinomial 
  distribution as the binomial and beta distributions are univariate versions 
  of the multinomial and Dirichlet distributions respectively. The special case 
  where α and β are integers is also known as the negative hypergeometric distribution. 

  ...

  To draw a beta-binomial random variate X ~ BetaBin(n,a,b) 
  simply draw p ~ Beta(a,b) and then draw X ~ B(n,p)
  """

  Below are some examples.

  Please also see ppl_distributions.pi and ppl_distributions_test.pi for
  more on this.

  Cf my Gamble model gamble_beta_binomial_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.
% import ordset.

main => go.


/*
  Here an example of beta_binomial(12,10,3)

  var : X
  Probabilities (truncated):
  11: 0.2438000000000000
  12: 0.2296000000000000
  10: 0.2040000000000000
  9: 0.1403000000000000
  .........
  4: 0.0050000000000000
  3: 0.0013000000000000
  2: 0.0002000000000000
  1: 0.0001000000000000
  mean = 10.0351


  theoretical_mean = 10.0


*/
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,
  println(theoretical_mean=beta_binomial_dist_mean(12,10,2)),


  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  N = 12,
  A = 10,
  B = 2,
  X = beta_binomial_dist(N,A,B),
  
  add("X",X).


/*

  From https://www.acsu.buffalo.edu/~adamcunn/probability/betabinomial.html
  """
  The probability each year that an apple contains a worm is a beta(0.5, 8) 
  random variable. 100 apples are picked one year. Let X be the number 
  containing worms.
  """

  var : X
  Probabilities (truncated):
  0: 0.2715000000000000
  1: 0.1237000000000000
  2: 0.0829000000000000
  3: 0.0708000000000000
  .........
  55: 0.0001000000000000
  53: 0.0001000000000000
  51: 0.0001000000000000
  50: 0.0001000000000000
  mean = 5.9199

  theoretical_mean = 5.88235
  

*/
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,
  println(theoretical_mean=beta_binomial_dist_mean(100,0.5,8)),


  % show_store_lengths,nl,
  % fail,
  nl.
go2 => true.
  
  
model2() =>
  A = 0.5,
  B = 8,
  X = beta_binomial_dist(100,A,B),
  
  add("X",X).


/*
  From https://www.acsu.buffalo.edu/~adamcunn/probability/betabinomial.html  
  """
  The probability that a random student can answer an exam question 
  correctly has a beta(3, 2) distribution. Let X be the number of correct 
  answers on a 20 question exam.
  """

  var : X
  Probabilities (truncated):
  14: 0.0795000000000000
  13: 0.0789000000000000
  12: 0.0773000000000000
  15: 0.0768000000000000
  .........
  3: 0.0175000000000000
  2: 0.0107000000000000
  1: 0.0061000000000000
  0: 0.0017000000000000
  mean = 12.0222

  theoretical_mean = 12.0


*/
go3 ?=>
  reset_store,
  run_model(10_000,$model3,[show_probs_trunc,mean % ,
                           % show_percentiles,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  println(theoretical_mean=beta_binomial_dist_mean(20,3,2)),


  % show_store_lengths,nl,
  % fail,
  nl.
go3 => true.
  
  
model3() =>
  A = 3,
  B = 2,
  X = beta_binomial_dist(20,A,B),
  
  add("X",X).

/*

  From https://www.acsu.buffalo.edu/~adamcunn/probability/betabinomial.html  
  """
  The probability that an archer in the Night's Watch can hit a white 
  walker has a beta(10, 3) distribution. A randomly chosen archer fires 
  36 arrows. Let X be the number of white walkers hit.
  """

  var : X
  Probabilities (truncated):
  30: 0.0930000000000000
  27: 0.0833333333333333
  31: 0.0823333333333333
  29: 0.0800000000000000
  .........
  12: 0.0013333333333333
  11: 0.0006666666666667
  8: 0.0006666666666667
  10: 0.0003333333333333
  mean = 27.6907


  theoretical_mean = 27.6923

*/
go4 ?=>
  reset_store,
  run_model(3_000,$model4,[show_probs_trunc,mean % ,
                           % show_percentiles,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  println(theoretical_mean=beta_binomial_dist_mean(36,10,3)),


  % show_store_lengths,nl,
  % fail,
  nl.
go4 => true.
  
  
model4() =>
  N = 36,
  A = 10,
  B = 3,
  X = beta_binomial_dist(N,A,B),
  
  add("X",X).