/* 

  Expected breakup length in Picat.

  https://brainstellar.com/puzzles/probability/204
  """
  A stick is broken into 3 pieces, by randomly choosing two points along 
  its unit length, and cutting it. What is the expected length of the middle part?

  Answer: 1/3
  """


  Cf my Gamble model gamble_expected_breakup_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.

main => go.

/*
  Using beta dist

  var : middle
  Probabilities (truncated):
  0.994182670195576: 0.0001000000000000
  0.991930528979104: 0.0001000000000000
  0.98845044973606: 0.0001000000000000
  0.987934127745989: 0.0001000000000000
  .........
  0.000329956244189: 0.0001000000000000
  0.000192992352875: 0.0001000000000000
  0.000174636291594: 0.0001000000000000
  0.000059289276112: 0.0001000000000000
  mean = 0.332011
  Percentiles:
  (0.001 0.00052593880289934661)
  (0.01 0.004669815253294388)
  (0.025 0.011634242473860648)
  (0.05 0.024645062352416994)
  (0.1 0.049460124065924878)
  (0.25 0.13101910205375022)
  (0.5 0.28827994896131393)
  (0.75 0.49844949201387456)
  (0.84 0.6013652152798814)
  (0.9 0.67972009053267413)
  (0.95 0.77736120617610616)
  (0.975 0.84504561152819957)
  (0.99 0.90215349763528807)
  (0.999 0.97278150480904368)
  (0.9999 0.99193075419322374)
  (0.99999 0.99395747859534334)


*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean,show_percentiles]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  P1 = beta_dist(1,1),
  P2 = beta_dist(1,1),
  Middle = abs(P1-P2),
  
  add("middle",Middle).


/*
  Using Dirichlet

  var : middle
  Probabilities (truncated):
  0.378536273321561: 0.0001000000000000
  0.377919683337133: 0.0001000000000000
  0.375122709917989: 0.0001000000000000
  0.37222125981185: 0.0001000000000000
  .........
  0.289254093696959: 0.0001000000000000
  0.28918048901601: 0.0001000000000000
  0.288992996356077: 0.0001000000000000
  0.283527445272028: 0.0001000000000000
  mean = 0.332449
  Percentiles:
  (0.001 0.29236407972297451)
  (0.01 0.30336630257494041)
  (0.025 0.30882108271519904)
  (0.05 0.31300062368888304)
  (0.1 0.31764968134600102)
  (0.25 0.3255834026336652)
  (0.5 0.33278070146988226)
  (0.75 0.33974194546738407)
  (0.84 0.34314011741598122)
  (0.9 0.34658727923073657)
  (0.95 0.35085617410240294)
  (0.975 0.35501961262968429)
  (0.99 0.35943971962053484)
  (0.999 0.36977047340247488)
  (0.9999 0.37791974499613068)
  (0.99999 0.37847462048901914)


*/
go2 ?=>
  reset_store,
  run_model(10_000,$model2,[show_probs_trunc,mean,show_percentiles]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go2 => true.
  
  
model2() =>
  % Note: We have to tweak the values of alpha.
  % An alpha of (vector 1 1 1) gives a too low value ~ 0.27
  Stick = dirichlet_dist([100,100,100]).sort,

  Middle = Stick[2],
  
  add("middle",Middle).