/* 

  Binomial process in Picat.

  From Mathematica BinomialProcess

  For more, see ppl_distributions.pi and ppl_distributions_test.pi


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


/*
  [p = 0.2,t = 10]
  PDF:
   2 0.301989888 
   1 0.268435456 
   3 0.201326592 
   0 0.1073741824 
   4 0.088080384 
   5 0.0264241152 
   6 0.005505024 
   7 0.000786432 
   8 0.000073728 
   9 0.000004096 

  CDF:
   0 0.1073741824 
   1 0.268435456 
   2 0.301989888 
   3 0.201326592 
   4 0.088080384 
   5 0.0264241152 
   6 0.005505024 
   7 0.000786432 
   8 0.000073728 
   9 0.000004096 

   Quantile:
   0.000001  0 
   0.00001  0 
   0.001  0 
   0.01  0 
   0.025  0 
   0.05  0 
   0.1  0 
   0.25  1 
   0.5  2 
   0.75  3 
   0.84  3 
   0.9  4 
   0.975  5 
   0.99  5 
   0.999  6 
   0.99999  8 
   0.999999  9 

   mean = 2.0
   variance = 1.6

*/
go ?=>
  P = 0.2,
  T = 10,
  println([p=P,t=T]),
  println("PDF:"),
  pdf_all($binomial_process_dist(P,T),0.01,0.999999).sort_down(2).printf_list,
  nl,
  println("CDF:"),
  pdf_all($binomial_process_dist(P,T),0.01,0.999999).printf_list,
  nl,
  println("Quantile:"),
  quantile_all($binomial_process_dist(P,T)).printf_list,
  nl,
  println(mean=binomial_process_dist_mean(P,T)),
  println(variance=binomial_process_dist_variance(P,T)),  
  nl.
go => true.
  
  
/*
  From Mathematica BinomialProcess
  """
  A quality assurance inspector randomly selects a series of 10 parts from a 
  manufacturing process that is known to produce 20% bad parts. Find the probability 
  that the inspector gets at most one bad part:

  selectionProcess = BinomialProcess[0.2];
  NProbability[x[10] <= 1, x \[Distributed] selectionProcess]
  -> 
  0.37581
  """

  Picat> X = binomial_process_dist_cdf(0.2,10,1)
  X = 0.375809638400001

  Picat> pdf_all($binomial_process_dist(0.2,10),0.01,0.999999).sort_down(2).printf_list
   2 0.301989888 
   1 0.268435456 
   3 0.201326592 
   0 0.1073741824 
   4 0.088080384 
   5 0.0264241152 
   6 0.005505024 
   7 0.000786432 
   8 0.000073728 
   9 0.000004096 


  var : val
  Probabilities:
  2: 0.3050000000000000
  1: 0.2653000000000000
  3: 0.2011000000000000
  0: 0.1135000000000000
  4: 0.0815000000000000
  5: 0.0258000000000000
  6: 0.0064000000000000
  7: 0.0013000000000000
  8: 0.0001000000000000
  mean = 1.9819

  var : prob
  Probabilities:
  false: 0.6212000000000000
  true: 0.3788000000000000
  mean = 0.3788

  theoretical_mean = 2.0
  theoretical_prob = 0.37581



*/
go2 ?=>
  P = 2/10,
  T = 10,
  Check = 1,
  reset_store,
  run_model(10_000,$model2(P,T,Check),[show_probs_trunc,mean % ,
                           % show_simple_stats % ,
                           % show_percentiles,
                           % show_hpd_intervals,hpd_intervals=[0.84,0.9,0.94,0.99,0.99999],
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  println(theoretical_mean=binomial_process_dist_mean(P,T)),                          
  println(theoretical_prob=binomial_process_dist_cdf(P,T,Check)),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go2 => true.
  
  
model2(P,T,Check) =>
  X = binomial_dist_n(T,P, T),

  Val = X.last,
  Prob = check(Val <= Check),
  add("val",Val),
  add("prob",Prob).