/* 

  Coin competition in Picat.

  From Pascal Bercker: "Coin Competition and the Geometric Distribution"
  https://medium.com/@pbercker/coin-competition-and-the-geometric-distribution-0f68ef440ea8 
  """
  Ty and Guy are both flipping fair coins until they respectively obtain 
  their first heads. Find the probability that it takes Ty at least 4 times as many 
  flips to obtain his first heads as Guy.
  """

  Cf my Gamble model gamble_coin_competition.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 : ty
  Probabilities:
  1: 0.4982500000000000
  2: 0.2508000000000000
  3: 0.1242000000000000
  4: 0.0641500000000000
  5: 0.0310500000000000
  6: 0.0157500000000000
  7: 0.0078000000000000
  8: 0.0037000000000000
  9: 0.0021000000000000
  10: 0.0010500000000000
  11: 0.0006500000000000
  13: 0.0002000000000000
  12: 0.0001500000000000
  14: 0.0001000000000000
  15: 0.0000500000000000
  mean = 2.0061

  var : guy
  Probabilities:
  1: 0.4980500000000000
  2: 0.2509500000000000
  3: 0.1253500000000000
  4: 0.0630500000000000
  5: 0.0332500000000000
  6: 0.0146500000000000
  7: 0.0076000000000000
  8: 0.0033000000000000
  9: 0.0017500000000000
  10: 0.0011500000000000
  11: 0.0005500000000000
  14: 0.0001500000000000
  12: 0.0001500000000000
  13: 0.0000500000000000
  mean = 1.9998

  var : p
  Probabilities:
  false: 0.9340000000000001
  true: 0.0660000000000000
  mean = 0.066

*/
go ?=>
  reset_store,
  run_model(20_000,$model,[show_probs,mean]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  

f(N,C,Target) = Res =>
  if C == Target then
    Res = N
  else
    T = bern(1/2),
    Res = f(N+1,C+T,Target)
  end.
  
model() =>
  Target = 1,

  Ty = f(0,0,Target),
  Guy = f(0,0,Target),

  % probability that it takes Ty at least 4 times as many flips to obtain his first heads as Guy.
  P = check(Ty >= 4*Guy),
  
  add("ty",Ty),
  add("guy",Guy),
  add("p",P).