/* 

  Rolling dice problem in Picat.

  From https://dtai.cs.kuleuven.be/problog/tutorial/basic/03_dice.html
  """
  Let’s consider an infinite number of dice, which we roll one after the other until we 
  see a six for the first time. What is the probability of stopping after n dice? 
  The first die is always rolled, those with higher numbers D are only rolled if the 
  previous roll did not stop the process.
  """

  Cf my Gamble model gamble_rolling_dice_problem4.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.

% Closed form of the probability that the process stops
% after N rolls.
theoretical(N) = (1/6) * (5/6)**(N-1).

/*
  Simulation:

  var : t
  Probabilities:
  1: 0.1679000000000000
  2: 0.1367000000000000
  3: 0.1164000000000000
  4: 0.0982000000000000
  5: 0.0786000000000000
  6: 0.0682000000000000
  7: 0.0542000000000000
  8: 0.0463000000000000
  9: 0.0383000000000000
  10: 0.0331000000000000
  11: 0.0261000000000000
  12: 0.0224000000000000
  13: 0.0174000000000000
  14: 0.0154000000000000
  15: 0.0130000000000000
  16: 0.0119000000000000
  17: 0.0106000000000000
  18: 0.0085000000000000
  20: 0.0055000000000000
  19: 0.0053000000000000
  21: 0.0045000000000000
  22: 0.0037000000000000
  24: 0.0024000000000000
  25: 0.0021000000000000
  23: 0.0020000000000000
  26: 0.0017000000000000
  28: 0.0016000000000000
  27: 0.0015000000000000
  30: 0.0008000000000000
  29: 0.0008000000000000
  33: 0.0007000000000000
  32: 0.0007000000000000
  35: 0.0006000000000000
  31: 0.0006000000000000
  34: 0.0005000000000000
  36: 0.0003000000000000
  47: 0.0002000000000000
  40: 0.0002000000000000
  39: 0.0002000000000000
  38: 0.0002000000000000
  37: 0.0002000000000000
  51: 0.0001000000000000
  49: 0.0001000000000000
  44: 0.0001000000000000
  42: 0.0001000000000000
  41: 0.0001000000000000
  mean = 6.0225

  The closed form (see theoretical/1 above):

  Picat> foreach(V in 1..20) println(V=((1/6) * (5/6)**(V-1))) end
  1 = 0.166667
  2 = 0.138889
  3 = 0.115741
  4 = 0.0964506
  5 = 0.0803755
  6 = 0.0669796
  7 = 0.0558163
  8 = 0.0465136
  9 = 0.0387613
  10 = 0.0323011
  11 = 0.0269176
  12 = 0.0224313
  13 = 0.0186928
  14 = 0.0155773
  15 = 0.0129811
  16 = 0.0108176
  17 = 0.00901465
  18 = 0.00751221
  19 = 0.00626017
  20 = 0.00521681


  Cf with the geometric distribution: geometric_dist(1/6) which 
  - in PPL - is the number of _failures_ before success:

  Picat> foreach(V in 0..20) println(V=geometric_dist_pdf(1/6,V)) end
  0 = 0.166667
  1 = 0.138889
  2 = 0.115741
  3 = 0.0964506
  4 = 0.0803755
  5 = 0.0669796
  6 = 0.0558163
  7 = 0.0465136
  8 = 0.0387613
  9 = 0.0323011
  10 = 0.0269176
  11 = 0.0224313
  12 = 0.0186928
  13 = 0.0155773
  14 = 0.0129811
  15 = 0.0108176
  16 = 0.00901465
  17 = 0.00751221
  18 = 0.00626017
  19 = 0.00521681
  20 = 0.00434734

*/
go ?=>
  N = 6,
  reset_store,
  run_model(10_000,$model(N),[show_probs,mean]),
  nl,
  % show_store_lengths,
  % fail,
  nl.
go => true.
  
roll(T,N) = Ret =>
  R = random_integer1(N),
  if R == N then
    Ret = T
  else
    Ret = roll(T+1,N)
  end.
% Another approach
roll2(T,N) = cond(R==N,T,roll2(T+1,N)) => R = random_integer1(N).  
model(N) =>
  T = roll2(1,N),
  add("t",T).