/* 

  Frustration patience in Picat.

  Grinstead & Snells "Introduction of Probability", page 86
  (https://math.dartmouth.edu/~prob/prob/prob.pdf)
  """
  Recently, a related problem appeared in a column of Marilyn vos Savant.
  Charles Price wrote to ask about his experience playing a certain form of solitaire,
  sometimes called 'frustration solitaire'. In this particular game, a deck of cards
  is shuffled, and then dealt out, one card at a time. As the cards are being dealt,
  the player counts from 1 to 13, and then starts again at 1. (Thus, each number is
  counted four times.) If a number that is being counted coincides with the rank of
  the card that is being turned up, then the player loses the game. Price found that
  he rarely won and wondered how often he should win. Vos Savant remarked that
  the expected number of matches is 4 so it should be difficult to win the game.
  """

  The expected number is 4, but it's quite a range of possible number of matches:

  Cf my Gamble model gamble_frustration_patience.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 : s
  Probabilities:
  3: 0.2036000000000000
  4: 0.2010000000000000
  5: 0.1631000000000000
  2: 0.1423000000000000
  6: 0.1047000000000000
  1: 0.0717000000000000
  7: 0.0552000000000000
  8: 0.0252000000000000
  0: 0.0159000000000000
  9: 0.0101000000000000
  10: 0.0050000000000000
  11: 0.0013000000000000
  12: 0.0005000000000000
  13: 0.0003000000000000
  14: 0.0001000000000000
  mean = 3.9693
  Percentiles:
  (0.001 0)
  (0.01 0)
  (0.025 1)
  (0.05 1)
  (0.1 2)
  (0.25 3)
  (0.5 4)
  (0.75 5)
  (0.84 6)
  (0.9 6)
  (0.95 7)
  (0.975 8)
  (0.99 9)
  (0.999 11)
  (0.9999 13.000099999999293)
  (0.99999 13.900010000001203)

  What is the probability of a win, i.e. no match:

  var : p
  Probabilities:
  false: 0.9841000000000000
  true: 0.0159000000000000
  mean = [false = 0.9841,true = 0.0159]
  Percentiles:
  (0.001 false)
  (0.01 false)
  (0.025 false)
  (0.05 false)
  (0.1 false)
  (0.25 false)
  (0.5 false)
  (0.75 false)
  (0.84 false)
  (0.9 false)
  (0.95 false)
  (0.975 false)
  (0.99 true)
  (0.999 true)
  (0.9999 true)
  (0.99999 true)

  
  The theoretical value is 1/exp(1)**4 = 0.01831564

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs,mean,
                           show_percentiles %, show_histogram,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  N = 52,
  A = 1..N,

  Cards = draw_without_replacement(N,A),
  Counts = [ cond(Cards[I] mod 13 == I mod 13,1,0)  : I in 1..N],

  S = Counts.sum,

  % What is the probability that I win, i.e. no match
  P = check(S == 0),
  
  add("s",S),
  add("p",P).