/* 

  Indistinguishable dice in Picat.

  From MIT 6.262 Discrete Stochastic Processes
  2. More Review; The Bernoulli Process
  https://www.youtube.com/watch?v=d4xfax4_Iww&list=PLEEF5322B331C1B98&index=2

  """
  Trivial example: Roll a white die and a red die.
  There are 36 sample outcomes, (i, j), 1 ≤ i, j ≤ 6,
  taken as equiprobable by symmetry.
  
  Roll 2 indistinguishable white dice. The white and
  red outcomes (i, j) and (j, i) for i =
  j are now indistinguishable. There are now 21 ‘finest grain’
  outcomes, but no sane person would use these as
  sample points.
  """

  So we have two indistinguishable dice where (i,j) = (j,i). 
  Which we represents as a condition that i<=j.
  There are 21 possible outcomes.
  What are the probabilities of the sums?


  Cf my Gamble model gamble_indistinguishable_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 : d1+d2 (normal dice)
  Probabilities:
  7: 0.1642500000000000
  8: 0.1416000000000000
  6: 0.1394000000000000
  5: 0.1115000000000000
  9: 0.1084000000000000
  10: 0.0837500000000000
  4: 0.0817500000000000
  11: 0.0580500000000000
  3: 0.0573000000000000
  2: 0.0274000000000000
  12: 0.0266000000000000
  mean = 7.001

  var : d1 + d2 with d1 <= d2
  Probabilities:
  6: 0.1444730077120823
  8: 0.1434447300771208
  7: 0.1404455869751499
  5: 0.0982005141388175
  9: 0.0939160239931448
  10: 0.0929734361610968
  4: 0.0928877463581834
  11: 0.0518423307626392
  3: 0.0492716366752356
  2: 0.0469580119965724
  12: 0.0455869751499572
  mean = 6.99409

*/
go ?=>
  reset_store,
  run_model(20_000,$model,[show_probs,mean % ,
                           % show_percentiles,show_histogram,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  D1 = dice(),
  D2 = dice(),

  observe(D1 <= D2),

  add("d1+d2 (normal dice)",(D1 + D2)),

  if observed_ok then
    add("d1 + d2 with d1 <= d2",(D1+D2)),
  end.