/* 

  Six districts, six robberies in Picat.

  From https://medium.com/puzzle-sphere/six-districts-six-robberies-can-you-solve-this-problem-92d851696da1
  """
  A city with 6 districts has 6 robberies in a particular week. Assume the 
  robberies are located randomly, with all possibilities for which robbery 
  occurred where equally likely. What is the probability that some district 
  had more than 1 robbery?

  ...

  From now on, the solution depends on how we consider the robberies. If 
  they are 'indistinguishable', meaning that we only know that a robbery 
  occurred but don’t know when and how — like assigning identical stickers 
  - the event d'“111111' is unique and we can use the Bose-Einstein 
  formula to count the total number of combinations:

              6!5!
  P(not R) :  ---
              11!

  -> P(R) = 1- P(not R) ... ~ 98.78%


  If the robberies are 'distinguishable', for example by their date or other 
  means, the event '111111' is not unique because we can permute each robbery 
  across the 6 districts. Therefore, we also need to count all possible 
  distinguishable distributions of 6 robberies across 6 districts, yielding
                 6!
     P(not R) = ---
                6^6

  -> P(R) = 1 - P(not R) ... ~ 98.46%

  """


  Cf my Gamble model gamble_six_districts_six_robberies.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:
  true: 0.9845000000000000
  false: 0.0155000000000000
  mean = 0.9845


*/
go ?=>
  reset_store,
  run_model(10_000,$model,[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
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  N = 6,
  % In what district does robberies[i] occurred?
  Robberies = random_integer1_n(N,N),

  Districts = collect(Robberies),

  % Probablity that (at least) one district has more than one robbery
  P = check( Districts.keys.len < N ),
  
  add("p",P).
  % add("robberies",Robberies),
  % add("districts",Districts).