/* 

  Negative Binomial basketball in Picat.

  From Mathematica NegativeBinomialDistribution
  """ 
  A basketball player shoots free throws until he hits 4 of them. His probability 
  of scoring in any one of them is 0.7. 

  Simulate the process:
  ...

  Find the number of shots the player is expected to shoot:
  > Expectation[x + 4, x => NegativeBinomialDistribution[4, 0.7]]
  5.71429

  Find the probability that the player requires 4 shots:
  > NProbability[x + 4 == 4, x => NegativeBinomialDistribution[4, 0.7]]
  0.2401
  """

  Simulate the process:
  Picat> X=negative_binomial_dist_n(4,0.7,20)  
  X = [1,1,0,0,0,1,1,1,3,1,2,0,2,2,2,4,2,0,4,3]

  Find the number of shots the player is expected to shoot:
  Picat> X=4+negative_binomial_dist_mean(4,0.7)
  X = 5.714285714285714
  
  Find the probability that the player requires 4 shots:
  Picat> X=negative_binomial_dist_pdf(4,0.7, 0)
  X = 0.2401


  Cf my Gamble model gamble_negative_binomial_basketball.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.7572000000000000
  true: 0.2428000000000000
  mean = [false = 0.7572,true = 0.2428]

  var : x
  Probabilities (truncated):
  1: 0.2908000000000000
  0: 0.2428000000000000
  2: 0.2226000000000000
  3: 0.1251000000000000
  .........
  8: 0.0024000000000000
  9: 0.0015000000000000
  10: 0.0003000000000000
  11: 0.0002000000000000
  mean = 1.6804

  var : y
  Probabilities (truncated):
  5: 0.2908000000000000
  4: 0.2428000000000000
  6: 0.2226000000000000
  7: 0.1251000000000000
  .........
  12: 0.0024000000000000
  13: 0.0015000000000000
  14: 0.0003000000000000
  15: 0.0002000000000000
  mean = 5.6804

*/
go ?=>
  println("Simulate the proccess"=negative_binomial_dist_n(4,0.7,20)),
  
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean]),
  nl,
  % show_store_lengths,
  % fail,
  nl.
go => true.
  
  
model() =>
  BN = 4,
  BP = 0.7,

  X = negative_binomial_dist(BN,BP),

  % Number of actual throws needed
  Y = X + 4,

  % What is the probability that it requires 4 shots
  P = check(Y==4),
  
  add("x",X),
  add("y",Y),  
  add("p",P).