/* 

  The blind archer in Picat.

  https://brainstellar.com/puzzles/probability/1006
  """
  A very sharp, consistently skillful blind archer aimed for the 
  center of a circular board and shot 2 arrows. He is expected to 
  hit the aim, but doesn't hit it for sure. The archer is told that 
  his first shot is better than second. He tried one more shot. What is 
  the probability that this 3rd shot is the best shot among 3? (ie, 
  Probability that 3rd arrow lands closer to center than his first two shots?)

  Answer: 1/3

  Solution: ... Notice that if there were (N-1) arrows, with first being 
            better than rest, and then he shoots Nth arrow. The probability 
            that Nth shot is best is 1/N.
  """

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

/*
  * With precision (Sigma) = 1

  var : arrow1
  Probabilities (truncated):
  2.453541712321077: 0.0001991238550378
  2.310196196133302: 0.0001991238550378
  2.169009915032114: 0.0001991238550378
  1.989660932008389: 0.0001991238550378
  .........
  -2.092763091419937: 0.0001991238550378
  -2.159588249830571: 0.0001991238550378
  -2.283174141538086: 0.0001991238550378
  -2.31981364765812: 0.0001991238550378
  mean = -0.0104582
  HPD intervals:
  HPD interval (0.84): -0.88072545365920..0.78758195682231

  var : arrow2
  Probabilities (truncated):
  3.714530845084849: 0.0001991238550378
  3.459532635940997: 0.0001991238550378
  3.385931262108365: 0.0001991238550378
  3.248642312391634: 0.0001991238550378
  .........
  -3.341012583632112: 0.0001991238550378
  -3.681310265188174: 0.0001991238550378
  -3.768376202798677: 0.0001991238550378
  -4.489137832857447: 0.0001991238550378
  mean = 0.0272834
  HPD intervals:
  HPD interval (0.84): -1.66313772501039..1.78124163356558

  var : arrow3
  Probabilities (truncated):
  3.238176949708494: 0.0001991238550378
  3.179551218732803: 0.0001991238550378
  3.072169595260815: 0.0001991238550378
  3.072129722871604: 0.0001991238550378
  .........
  -3.247955288384437: 0.0001991238550378
  -3.283151705185248: 0.0001991238550378
  -3.290174281044576: 0.0001991238550378
  -3.509903774157841: 0.0001991238550378
  mean = -0.00652723
  HPD intervals:
  HPD interval (0.84): -1.39827316638491..1.41852836765855

  var : p
  Probabilities:
  false: 0.6652727996814018
  true: 0.3347272003185982
  mean = [false = 0.665273,true = 0.334727]
  HPD intervals:
  show_hpd_intervals: data is not numeric

  * With better precision (Sigma) = 0.1

  var : arrow1
  Probabilities (truncated):
  0.306823563355117: 0.0002003606491685
  0.253212427484149: 0.0002003606491685
  0.241758513176084: 0.0002003606491685
  0.228480073061236: 0.0002003606491685
  .........
  -0.207758296998262: 0.0002003606491685
  -0.209587565499239: 0.0002003606491685
  -0.225742251476341: 0.0002003606491685
  -0.262506764172004: 0.0002003606491685
  mean = 4.3035e-05
  HPD intervals:
  HPD interval (0.84): -0.08270120891570..0.08444643862243

  var : arrow2
  Probabilities (truncated):
  0.359106655594433: 0.0002003606491685
  0.352590578439617: 0.0002003606491685
  0.348092574827394: 0.0002003606491685
  0.337621794905395: 0.0002003606491685
  .........
  -0.338287047212688: 0.0002003606491685
  -0.342106399108857: 0.0002003606491685
  -0.355433075652612: 0.0002003606491685
  -0.364886474499959: 0.0002003606491685
  mean = 0.00185081
  HPD intervals:
  HPD interval (0.84): -0.16804335029013..0.17275907640797

  var : arrow3
  Probabilities (truncated):
  0.414599744763812: 0.0002003606491685
  0.379174847266274: 0.0002003606491685
  0.348142440379716: 0.0002003606491685
  0.345562377493318: 0.0002003606491685
  .........
  -0.326243468026795: 0.0002003606491685
  -0.326757985969109: 0.0002003606491685
  -0.343474152397267: 0.0002003606491685
  -0.354307232145705: 0.0002003606491685
  mean = -0.000588472
  HPD intervals:
  HPD interval (0.84): -0.13951592548208..0.14797319221927

  var : p
  Probabilities:
  false: 0.6754157483470247
  true: 0.3245842516529753
  mean = [false = 0.675416,true = 0.324584]
  HPD intervals:
  show_hpd_intervals: data is not numeric

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean,show_hpd_intervals,hpd_intervals=[0.84]]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  % Sigma = 1,
  Sigma = 0.1, % better precision
  Arrow1 = normal_dist(0,Sigma),
  Arrow2 = normal_dist(0,Sigma),
  Arrow3 = normal_dist(0,Sigma),  
  
  observe( abs(Arrow1) < abs(Arrow2)),

  P = check(abs(Arrow3) < abs(Arrow1)),
  
  if observed_ok then
    add("arrow1",Arrow1),
    add("arrow2",Arrow2),
    add("arrow3",Arrow3),
    add("p",P)
  end.


/*
  Different number of arrows

  num_arrows = 3
  var : best
  Probabilities:
  3: 0.3758907363420427
  2: 0.3640142517814727
  1: 0.2600950118764845
  mean = 2.1158
  HPD intervals:
  HPD interval (0.84): 1.00000000000000..3.00000000000000

  var : p
  Probabilities:
  false: 0.6241092636579573
  true: 0.3758907363420427
  mean = [false = 0.624109,true = 0.375891]
  HPD intervals:
  show_hpd_intervals: data is not numeric


  num_arrows = 7
  var : best
  Probabilities:
  6: 0.1775510204081633
  4: 0.1707482993197279
  2: 0.1639455782312925
  5: 0.1619047619047619
  3: 0.1612244897959184
  7: 0.1551020408163265
  1: 0.0095238095238095
  mean = 4.46463
  HPD intervals:
  HPD interval (0.84): 1.00000000000000..6.00000000000000

  var : p
  Probabilities:
  false: 0.8448979591836735
  true: 0.1551020408163265
  mean = [false = 0.844898,true = 0.155102]
  HPD intervals:
  show_hpd_intervals: data is not numeric


  num_arrows = 15
  var : best
  Probabilities (truncated):
  5: 0.0881889763779528
  9: 0.0803149606299213
  14: 0.0787401574803150
  15: 0.0771653543307087
  .........
  4: 0.0645669291338583
  13: 0.0598425196850394
  7: 0.0582677165354331
  10: 0.0551181102362205
  mean = 8.47874
  HPD intervals:
  HPD interval (0.84): 2.00000000000000..13.00000000000000

  var : p
  Probabilities:
  false: 0.9228346456692913
  true: 0.0771653543307087
  mean = [false = 0.922835,true = 0.0771654]
  HPD intervals:
  show_hpd_intervals: data is not numeric

*/ 
go2 ?=>
  member(NumArrows,[3,7,15]),
  println(num_arrows=NumArrows),
  reset_store,
  run_model(10_000,$model2(NumArrows),[show_probs_trunc,mean,show_hpd_intervals,hpd_intervals=[0.84]]),
  nl,
  % show_store_lengths,nl,
  fail,
  nl.
go2 => true.
  
  
model2(NumArrows) =>
  Sigma = 1,
  Arrows = normal_dist_n(0,Sigma,NumArrows),

  foreach(I in 2..NumArrows)
    observe( abs(Arrows[1]) < abs(Arrows[I]))
  end,

  Best = argmin(Arrows).first,

  P = check(Best == NumArrows),
  
  if observed_ok then
    add("best",Best),
    add("p",P)
  end.