/* 

  How long until you roll two consecutive sizes in Picat.

  Dr. Robert Kübler "How Long Until You Roll Two Consecutive Sixes?"
  https://www.cantorsparadise.com/how-long-until-you-roll-two-consecutive-sixes-e184292141d5
  """
  You roll a fair standard die repeatedly, until you get two consecutive sixes. 
  What is the expected number of rolls before the game ends?
  """

  Cf my Gamble model gamble_how_long_until_two_conecutive_sixes.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 : len
  Probabilities (truncated):
  2: 0.0299000000000000
  3: 0.0237000000000000
  7: 0.0233000000000000
  5: 0.0222000000000000
  11: 0.0207000000000000
  8: 0.0203000000000000
  4: 0.0203000000000000
  15: 0.0197000000000000
  9: 0.0197000000000000
  6: 0.0197000000000000
  .........
  226: 0.0001000000000000
  222: 0.0001000000000000
  218: 0.0001000000000000
  211: 0.0001000000000000
  207: 0.0001000000000000
  202: 0.0001000000000000
  201: 0.0001000000000000
  199: 0.0001000000000000
  195: 0.0001000000000000
  194: 0.0001000000000000
  mean = 42.1928
  [len = 10000,min = 2,mean = 42.1928,median = 30.0,max = 373,variance = 1672.8,stdev = 40.8999]
  HPD intervals:
  HPD interval (0.84): 2.00000000000000..76.00000000000000
  HPD interval (0.9): 2.00000000000000..95.00000000000000
  HPD interval (0.94): 2.00000000000000..116.00000000000000
  HPD interval (0.99): 2.00000000000000..192.00000000000000
  HPD interval (0.99999): 2.00000000000000..373.00000000000000


*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean,truncate_size=10,
                           show_simple_stats,
                           % show_percentiles,
                           show_hpd_intervals,hpd_intervals=[0.84,0.9,0.94,0.99,0.99999]
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
f(A) = Res =>
  T = A.take_last(2),
  if T == [6,6] then
    Res = A
  else
    Res = f(A++[dice()])
  end.
  
model() =>
  A = f([]),
  Len = A.len,
  
  add("len",Len).