/* 

  Newton-Pepy's problem in Picat.

  https://en.wikipedia.org/wiki/Newton%E2%80%93Pepys_problem
  """
  The Newton–Pepys problem is a probability problem concerning the probability of throwing 
  sixes from a certain number of dice.

  In 1693 Samuel Pepys and Isaac Newton corresponded over a problem posed to Pepys by a school 
  teacher named John Smith. The problem was:

  Which of the following three propositions has the greatest chance of success?

  A. Six fair dice are tossed independently and at least one "6" appears.
  B. Twelve fair dice are tossed independently and at least two "6"s appear.
  C. Eighteen fair dice are tossed independently and at least three "6"s appear.

  Pepys initially thought that outcome C had the highest probability, but Newton correctly concluded 
  that outcome A actually has the highest probability.

  Solution
  The probabilities of outcomes A, B and C are:

  P(A): 0.6651
  P(B): 0.6187
  P(C): 0.5973
  """

  Cf my Gamble model gamble_newton_pepys_problem.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 : a
  Probabilities:
  true: 0.6667999999999999
  false: 0.3332000000000000
  mean = 0.6668

  var : b
  Probabilities:
  true: 0.6201000000000000
  false: 0.3799000000000000
  mean = 0.6201

  var : c
  Probabilities:
  true: 0.6027000000000000
  false: 0.3973000000000000
  mean = 0.6027

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean % ,
                           % show_percentiles,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  A = check(binomial_dist(6,1/6)  >= 1),
  B = check(binomial_dist(12,1/6) >= 2),
  C = check(binomial_dist(18,1/6) >= 3),  
  add("a",A),
  add("b",B),
  add("c",C).



/*
  Basic simulation

  var : a
  Probabilities:
  true: 0.6596000000000000
  false: 0.3404000000000000
  mean = 0.6596

  var : b
  Probabilities:
  true: 0.6266000000000000
  false: 0.3734000000000000
  mean = 0.6266

  var : c
  Probabilities:
  true: 0.6058000000000000
  false: 0.3942000000000000
  mean = 0.6058

*/
go2 ?=>
  reset_store,
  run_model(10_000,$model2,[show_probs_trunc,mean % ,
                           % show_percentiles,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go2 => true.
  
  
model2() =>
  A = check( [ cond(dice() == 6,1,0) : _ in 1.. 6].sum  >= 1),
  B = check( [ cond(dice() == 6,1,0) : _ in 1..12].sum  >= 2),
  C = check( [ cond(dice() == 6,1,0) : _ in 1..18].sum  >= 3),
  add("a",A),
  add("b",B),
  add("c",C).