/* 

  Minimum (and maximum) of 4 dice in Picat.

  From Mathematica OrderDistribution
  """
  Four six-sided dice are rolled. Find the expectation of the minimum value:
  -->
  2275/1296
  1.7554

  Find the expectation of the maximum value:
  -->
  6797/1296
  5.2446
  """


  Cf my Gamble model gamble_minimum_of_four_dice.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 : min val
  Probabilities:
  1: 0.5250000000000000
  2: 0.2832000000000000
  3: 0.1305000000000000
  4: 0.0477000000000000
  5: 0.0127000000000000
  6: 0.0009000000000000
  mean = 1.7426

  var : max val
  Probabilities:
  6: 0.5179000000000000
  5: 0.2791000000000000
  4: 0.1381000000000000
  3: 0.0512000000000000
  2: 0.0127000000000000
  1: 0.0010000000000000
  mean = 5.2353


  order_statistics_with_replacement_discrete_pdf(PDF,CDF,N,1,I)
  Min value:
  1 = 0.517747
  2 = 0.284722
  3 = 0.135031
  4 = 0.0501543
  5 = 0.0115741
  6 = 0.000771605

  Max value:
  1 = 0.000771605
  2 = 0.0115741
  3 = 0.0501543
  4 = 0.135031
  5 = 0.284722
  6 = 0.517747


*/
go ?=>
  N = 4,
  reset_store,
  run_model(10_000,$model(N),[show_probs_trunc,mean % ,
                           % show_percentiles,show_histogram,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  println("order_statistics_with_replacement_discrete_pdf(PDF,CDF,N,1,I):"),
  PDF=$random_integer1_pdf(6),
  CDF=$random_integer1_cdf(6),
  println("Min value:"),
  foreach(I in 1..6) 
    println(I=order_statistics_with_replacement_discrete_pdf(PDF,CDF,N,1,I))
  end,
  println("\nMax value:"),
  foreach(I in 1..6) 
    println(I=order_statistics_with_replacement_discrete_pdf(PDF,CDF,N,4,I))
  end,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model(N) =>
  Rolls = random_integer1_n(6,N),
  MinVal = Rolls.min,
  MaxVal = Rolls.max,  
  
  add("min val",MinVal),
  add("max val",MaxVal).