/* 

  Multinomial callcenter in Picat.

  From Mathematica (MultinomialDistribution)
  """
  In calling a customer service center, one of three things may happen: the line 
  is busy with probability 0.4, a caller gets the wrong party with probability 
  0.1, or a caller gets connected to an agent. Find the probability that a caller 
  calling at 6 different times gets a busy signal 4 times and twice connects 
  directly to an agent:

  D = MultinomialDistribution(6, (0.4, 0.1, 0.5))
  PDF(D, (4, 0, 2))
  -> 0.096
  """

  Cf my Gamble model gamble_multinomial_callcenter.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 : p
  Probabilities:
  false: 0.9035000000000000
  true: 0.0965000000000000
  mean = [false = 0.9035,true = 0.0965]

  var : p2
  Probabilities:
  false: 0.9069000000000000
  true: 0.0931000000000000
  mean = [false = 0.9069,true = 0.0931]

  Theoretical: 0.09600000000000


*/
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,
  printf("Theoretical: %.14f\n",multinomial_dist_pdf(6, [4/10,1/10,5/10],[4,0,2])),
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  N = 6,

  Calls = categorical_dist_n([4/10,1/10,5/10],[busy,wrong_party,connected],N),

  NumBusy = count_occurrences(busy,Calls),
  WrongParty = count_occurrences(wrong_party,Calls),
  Connected = count_occurrences(connected,Calls),    
  
  P = check([NumBusy,WrongParty,Connected] == [4,0,2]),
  P2 = check(multinomial_dist(6, [4/10,1/10,5/10])==[4,0,2]),

  add("p",P),
  add("p2",P2).