/* 

  Same rank in Picat.

  From Gunnar Blom, Lars Holst, Dennis Sandell:
  "Problems and Snapshots from the World of Probability"
  Page 22.
  """
  100 students are seating (randomly) on seats numbered 0 to 99. 
  The lectures says: assume that I rank you according to your age. 
  How many will have the same rank number as the seat number?
  """

  Below are some model showing the same expectation of same
  rank: E[same rank] = 1.
  Also, the probability of exactly 0 matching is about the same
  as the probability of (exactly) 1 match: about 0.3678 (1/exp(1)).
  

  Cf 
  - ppl_matching_distribution.pi
  - my Gamble model gamble_same_rank.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 : num same rank
  Probabilities:
  1: 0.3707000000000000
  0: 0.3620000000000000
  2: 0.1894000000000000
  3: 0.0605000000000000
  4: 0.0140000000000000
  5: 0.0028000000000000
  6: 0.0005000000000000
  7: 0.0001000000000000
  mean = 1.0047

  var : p
  Probabilities:
  true: 0.6380000000000000
  false: 0.3620000000000000
  mean = [true = 0.638,false = 0.362]

  The theoretical values:
  Picat> [V=matching_dist_pdf(100,V) : V in 0..20].print_list
  0 = 0.367879
  1 = 0.367879
  2 = 0.18394
  3 = 0.0613132
  4 = 0.0153283
  5 = 0.00306566
  6 = 0.000510944
  7 = 7.2992e-05
  8 = 9.12399e-06
  9 = 1.01378e-06
  10 = 1.01378e-07
  11 = 9.21616e-09
  12 = 7.68013e-10
  13 = 5.90779e-11
  14 = 4.21985e-12
  15 = 2.81323e-13
  16 = 1.75827e-14
  17 = 1.03428e-15
  18 = 0.0
  19 = 0.0
  20 = 0.0

  Picat> X=matching_dist_mean(100)                           
  X = 1

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean]),
  nl,
  % show_store_lengths,
  % fail,
  nl.
go => true.
  
  
model() =>
  N = 100,
  % X = shuffle(1..N),
  X = draw_without_replacement(N,1..N),
  NumSameRank = sum([1 : I in 1..N, X[I] == I]),

  P = check(NumSameRank > 0),
  
  add("num same rank",NumSameRank),
  add("p",P).