/* 

  Three men with hats in Picat.

  From Introduction to Probability Models by Sheldon Ross
  page 8
  """
  Example 1.8 
  Suppose that each of three men at a party throws his hat into the center 
  of the room. The hats are first mixed up and then each man randomly 
  selects a hat. What is the probability that none of the three men 
  selects his own hat?
  """

  Cf my Gamble model gamble_three_men_with_hats.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 : a
  Probabilities:
  2: 0.3535839160839161
  3: 0.3513986013986014
  1: 0.2950174825174825
  mean = 2.05638

  var : b
  Probabilities:
  1: 0.3461538461538461
  3: 0.3286713286713286
  2: 0.3251748251748252
  mean = 1.98252

  var : c
  Probabilities:
  1: 0.3588286713286714
  2: 0.3212412587412588
  3: 0.3199300699300699
  mean = 1.9611

  var : none gets their own hat
  Probabilities:
  false: 0.6472902097902098
  true: 0.3527097902097902
  mean = 0.35271

  var : at least one get their own hat
  Probabilities:
  true: 0.6472902097902098
  false: 0.3527097902097902
  mean = 0.64729

  var : exactly one gets their own hat
  Probabilities:
  true: 0.5008741258741258
  false: 0.4991258741258741
  mean = 0.500874

*/
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() =>
  Ps = [1/3,1/3,1/3],
  Vs = 1..3,

  L = categorical_n(Ps,Vs,3),
  [A,B,C] = L,
  
  observe(all_different(L)),

  NumGetOwnHat = [cond(L[I] == I,1,0) : I in 1..3].sum,

  P = check(NumGetOwnHat == 0), % None get their own hat: 2/6=1/3
  P2 = check(NumGetOwnHat >= 1), % At least one find their own hat: 2/3
  P3 = check(NumGetOwnHat == 1), % Exactly one get their own hat

  if observed_ok then
    add("a",A),
    add("b",B),
    add("c",C),
    add("none gets their own hat",P),
    add("at least one get their own hat",P2),
    add("exactly one gets their own hat",P3)
  end.