/* 

  Coin flips problem in Picat.

  From Cliff Pickover: https://x.com/pickover/status/1947785369253970055
  """
  Monica flips a fair coin until more head than tails have been flipped.
  What is the expected number of times of flipping the coin?
  """
  
  Cf my Gamble model gamble_coin_flips2.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.

/*

  p = 0.5
  var : total losses
  Probabilities:
  1: 0.5001000000000000
  2: 0.2481000000000000
  3: 0.1255000000000000
  4: 0.0635000000000000
  5: 0.0301000000000000
  6: 0.0177000000000000
  7: 0.0074000000000000
  8: 0.0036000000000000
  9: 0.0021000000000000
  11: 0.0011000000000000
  10: 0.0008000000000000
  mean = 2.0031

  var : num heads
  Probabilities:
  1: 0.5001000000000000
  2: 0.2481000000000000
  3: 0.1255000000000000
  4: 0.0635000000000000
  5: 0.0301000000000000
  6: 0.0177000000000000
  7: 0.0074000000000000
  8: 0.0036000000000000
  9: 0.0021000000000000
  10: 0.0014000000000000
  11: 0.0005000000000000
  mean = 2.0025

  var : num tails
  Probabilities:
  0: 0.5001000000000000
  1: 0.2481000000000000
  2: 0.1255000000000000
  3: 0.0635000000000000
  4: 0.0301000000000000
  5: 0.0177000000000000
  6: 0.0074000000000000
  7: 0.0036000000000000
  8: 0.0021000000000000
  9: 0.0008000000000000
  11: 0.0006000000000000
  10: 0.0005000000000000
  mean = 1.0037

  theoretical_mean = 2.0
  PDF:
   1 0.5 
   2 0.25 
   3 0.125 
   4 0.0625 
   5 0.03125 
   6 0.015625 
   7 0.0078125 
   8 0.00390625 
   9 0.001953125 
  10 0.0009765625 
  11 0.00048828125 
  12 0.000244140625 
  13 0.0001220703125 
  14 0.00006103515625 
  15 0.000030517578125 
  16 0.000015258789062 
  17 0.000007629394531 
  18 0.000003814697266 
  19 0.000001907348633 
  20 0.000000953674316 



*/
go ?=>
  P = 1/2,
  println(p=P),
  reset_store,
  run_model(10_000,$model(P),[show_probs,mean]),
  println(theoretical_mean=shifted_geometric_dist_mean(P)),                          
  println("PDF:"),
  [K=shifted_geometric_dist_pdf(P,K) : K in 1..20].sort_down(2).printf_list,
  
                          
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
f(N,H,T,P) = Res =>
  if N >= 11 ; H > T then
    Res = [N,H,T]
  else
    V = flip(1/2),
    N2 = N + cond(V==true,1,0),
    T2 = T + cond(V==false,1,0),    
    Res = f(N+1,N2,T2,P)
  end.
  
model(P) =>
  [TotalLosses,NumHeads,NumTails] = f(0,0,0,P),

  add("total losses",TotalLosses),
  add("num heads",NumHeads),
  add("num tails",NumTails).