/* 

  Negative binomal coins II in Picat.

  From Mathematica (NegativeBinomialDistribution)
  """
  A coin was flipped 10 times and the 8^th head occurred at 
  the 10^th flip. Find the probability of such an event if the coin is fair:

    fairD = NegativeBinomialDistribution(8, 1/2);
    Probability(x == 2, x fairD))
    -> 
      9/256 (0.0351563)

  Assuming the coin may not be fair, find the most likely value for p:
    D = NegativeBinomialDistribution(8, p);
    FindMaximum(PDF(D, 2), (p, 0.9))
    -> 
    (0.241592, (p -> 0.8))
  """

  * A coin was flipped 10 times and the 8^th head occurred at 
    the 10^th flip. Find the probability of such an event if the coin is fair:

    Picat> X=negative_binomial_dist_pdf(8,1/2,2)
    X = 0.03515625

  * Assuming the coin may not be fair, find the most likely value for p:
    Picat> X=[negative_binomial_dist_pdf(8,P,2) : P in 0..0.01..1],ArgMax=X.argmax.first,D=X[ArgMax]
    ...
    ArgMax = 81
    D = 0.2415919104

  Cf my Gamble model gamble_binomial_coins2.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 : c
  Probabilities:
  7: 0.1041000000000000
  6: 0.1039000000000000
  8: 0.0996000000000000
  5: 0.0926000000000000
  9: 0.0898000000000000
  4: 0.0781000000000000
  10: 0.0759000000000000
  11: 0.0599000000000000
  3: 0.0574000000000000
  12: 0.0506000000000000
  13: 0.0407000000000000
  2: 0.0335000000000000
  14: 0.0264000000000000
  15: 0.0192000000000000
  1: 0.0175000000000000
  16: 0.0161000000000000
  17: 0.0099000000000000
  18: 0.0078000000000000
  0: 0.0046000000000000
  19: 0.0041000000000000
  20: 0.0034000000000000
  21: 0.0017000000000000
  22: 0.0012000000000000
  23: 0.0009000000000000
  24: 0.0004000000000000
  27: 0.0002000000000000
  25: 0.0002000000000000
  30: 0.0001000000000000
  28: 0.0001000000000000
  26: 0.0001000000000000
  mean = 8.0244

  var : p
  Probabilities:
  false: 0.9665000000000000
  true: 0.0335000000000000
  mean = 0.0335

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs,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() =>
  C = negative_binomial_dist(8,1/2),
  P = check(C == 2),
  
  add("c",C),
  add("p",P).