/* 

  Negative binomial selling candies in Picat.

  From https://en.wikipedia.org/wiki/Negative_binomial_distribution
  """
  Selling candy
  Pat Collis is required to sell candy bars to raise money for the 6th grade field trip. 
  Pat is (somewhat harshly) not supposed to return home until five candy bars have been 
  sold. So the child goes door to door, selling candy bars. At each house, there is a 
  0.6 probability of selling one candy bar and a 0.4 probability of selling nothing.

  What's the probability of selling the last candy bar at the nth house?
  ...
  What's the probability that Pat finishes on the tenth house?
    f(10)=0.1003290624
  """

  Picat> X=negative_binomial_dist_mean(10-5,6/10)  
  X = 3.333333333333333

  * What's the probability of selling the last candy bar at the nth house?

    Picat> pdf_all($negative_binomial_dist(10-5,6/10),0.001,0.99999).printf_list
     0 0.07776 
     1 0.15552 
     2 0.186624 
     3 0.1741824 
     4 0.13934592 
     5 0.1003290624 
     6 0.0668860416 
     7 0.04204265472 
     8 0.025225592832 
     9 0.0145747869696 
    10 0.008161880702976 
    11 0.004451934928896 
    12 0.002374365295411 
    13 0.001241975692984 
    14 0.000638730356392 
    15 0.000323623380572 
    16 0.000161811690286 
    17 0.000079954011671 
    18 0.000039088627928 
    19 0.000018927125102 
    20 0.000009085020049 

  * What's the probability that Pat finishes on the tenth house?

    Picat> X=negative_binomial_dist_pdf(10-5,6/10,5)                                                
    X = 0.1003290624
  

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

/*

  var : y
  Probabilities (truncated):
  2: 0.1863000000000000
  3: 0.1722000000000000
  1: 0.1495000000000000
  4: 0.1412000000000000
  .........
  13: 0.0009000000000000
  14: 0.0006000000000000
  16: 0.0002000000000000
  15: 0.0002000000000000
  mean = 3.3336

  var : p
  Probabilities:
  false: 0.8927000000000000
  true: 0.1073000000000000
  mean = 0.1073

*/
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,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  NumHouses  = 10,
  NumCandies = 5,
  ProbOfSeeling = 6/10,

  Y = negative_binomial_dist(NumHouses-NumCandies,ProbOfSeeling),
  P = check(Y == NumCandies),
  
  add("y",Y),
  add("p",P).