/* 

  Rope length in Picat.

  From Wolfram U's Intro to Probability 
  https://www.wolframcloud.com/obj/online-courses/introduction-to-probability/joint-distributions.html
  """
  A rope of length 1 is cut according to a 15-10-5 Dirichlet distribution.
  ...
  Find the expected length of the rope parts:
  """

  Mathematica's solution:
  dist = DirichletDistribution[{15, 10, 5}]
  Table[Mean[MarginalDistribution[dist, i]], {i, 2}]
  {1/2, 1/3}


  Compare with
  Picat> X=dirichlet_dist_mean([15,10,5])
  X = [0.5,0.333333333333333,0.166666666666667]

  Picat> X=dirichlet_dist_mean([15,10,5]).map(to_rat)
  X = [1 / 2,1 / 3,1 / 6]

  Note that the numbers 15,10,5 represents the expected proportions:
  Picat> X=[15,10,5], P = [V/X.sum : V in X]
  P = [0.5,0.333333333333333,0.166666666666667]


  Cf my Gamble model gamble_rope_length.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 : rope 1
  mean = 0.500809

  var : rope 2
  mean = 0.332482

  var : rope 3
  mean = 0.166709

  var : len
  mean = 1.0

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[mean]),
  nl,
  nl.
go => true.
  
  
model() =>
  Rope = dirichlet_dist([15,10,5]),
  
  add("rope 1",Rope[1]),
  add("rope 2",Rope[2]),
  add("rope 3",Rope[3]),
  add("len",Rope.sum).

/*
  Compare with multinomial dist

  Picat> X=multinomial_dist_mean(100,[15,10,5].simplex)             
  X = [50.000000000000007,33.333333333333336,16.666666666666668]

  var : rope 1
  mean = 0.499552

  var : rope 2
  mean = 0.332466

  var : rope 3
  mean = 0.167982

  var : len
  mean = 1.0

*/
go2 ?=>
  reset_store,
  run_model(10_000,$model,[mean]),
  nl,
  nl.
go2 => true.
  
  
model2() =>
  Rope = multinomial_dist(100,[15,10,5].simplex),
  
  add("rope 1",Rope[1]),
  add("rope 2",Rope[2]),
  add("rope 3",Rope[3]),
  add("len",Rope.sum).