/* 

  Discrete uniform dist in Picat.

  From Mathematica DiscreteUniformDistribution

  See ppl_distributions.pi and ppl_distributions_test.pi for more on this.


  Cf my Gamble model gamble_discrete_uniform_dist.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 : d
  Probabilities:
  9: 0.1749000000000000
  7: 0.1696000000000000
  5: 0.1673000000000000
  10: 0.1662000000000000
  6: 0.1658000000000000
  8: 0.1562000000000000
  mean = 7.5042

  Theoretical:
  Picat> pdf_all($discrete_uniform_dist(5,10)).printf_list
   5 0.166666666666667 
   6 0.166666666666667 
   7 0.166666666666667 
   8 0.166666666666667 
   9 0.166666666666667 
  10 0.166666666666667 

  Picat> cdf_all($discrete_uniform_dist(5,10)).printf_list
   5 0.166666666666667 
   6 0.333333333333333 
   7 0.5 
   8 0.666666666666667 
   9 0.833333333333333 
  10 1.0 

  Picat> quantile_all($discrete_uniform_dist(5,10)).printf_list
   0.000001  5 
   0.00001  5 
   0.001  5 
   0.01  5 
   0.025  5 
   0.05  5 
   0.25  6 
   0.5  8 
   0.75  9 
   0.84 10 
   0.975 10 
   0.99 10 
   0.999 10 
   0.99999 10 
   0.999999 10 


*/
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 = 5,
  B = 10,
  D = discrete_uniform_dist(A,B),
  
  add("d",D).

/*
  From Mathematica DiscreteUniformDistribution
  """
  A fair six-sided die can be modeled using a DiscreteUniformDistribution: 

  dice = DiscreteUniformDistribution[{1, 6}]

  Generate 10 throws of the die:
  RandomVariate[dice, 10]
  -> {5, 2, 1, 5, 3, 6, 1, 5, 5, 2}

  Compute the probability that the sum of three dice values is less than 6:
  Probability[ x1 + x2 + x3 <= 6, {x1 e dice, x2 e dice, x3 e dice}]
  -> 5/54
  -> 0.0925926
  """

  Picat> X = discrete_uniform_dist_n(1,6,10)                   
  X = [2,6,4,5,5,3,1,4,4,6]

  Picat PPL has some short cuts for this: dice and random_integer1 (1..N)

  Picat> X = dice_n(6,10)                   
  X = [3,3,3,3,2,1,6,1,1,5]

  Picat> X = random_integer1_n(6,10)
  X = [2,2,6,5,3,4,1,4,4,5]

  var : s
  Probabilities:
  10: 0.1273000000000000
  11: 0.1235000000000000
  12: 0.1183000000000000
  9: 0.1176000000000000
  8: 0.0942000000000000
  13: 0.0909000000000000
  14: 0.0744000000000000
  7: 0.0699000000000000
  6: 0.0465000000000000
  15: 0.0463000000000000
  5: 0.0289000000000000
  16: 0.0260000000000000
  17: 0.0152000000000000
  4: 0.0134000000000000
  3: 0.0043000000000000
  18: 0.0033000000000000
  mean = 10.494

  var : p
  Probabilities:
  false: 0.9069000000000000
  true: 0.0931000000000000
  mean = [false = 0.9069,true = 0.0931]



*/
go2 ?=>
  reset_store,
  run_model(10_000,$model2,[show_probs,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() =>
  N = 3, % number of dice
  Ds = discrete_uniform_dist_n(1,6,N),
  S = Ds.sum,
  P = check(S <= 6),
  
  add("s",S),
  add("p",P).