/* 

  Doomsday in Picat.

  http://www.homoexcelsior.com/archive/transhuman/msg03904.html
  """
  Imagine that two big urns are put in front of you, and you know that 
  one of them contains ten balls and the other a million, but you are 
  ignorant as to which is which. You know the balls in each urn are 
  numbered 1, 2, 3, 4 ... etc. Now you take a ball at random from the 
  left urn, and it is number 7. Clearly, this is a strong indication 
  that that urn contains only ten balls. If originally the odds were 
  fifty-fifty, a swift application of Bayes' theorem gives you the 
  posterior probability that the left urn is the one with only ten 
  balls. (Pposterior (L=10) = 0.999990).

  But now consider the case where instead of the urns you have two 
  possible human races, and instead of balls you have individuals, 
  ranked according to birth order. As a matter of fact, you happen to 
  find that your rank is about sixty billion. Now, say Carter and Leslie, 
  we should reason in the same way as we did with the urns. That you 
  should have a rank of sixty billion or so is much more likely if 
  only 100 billion persons will ever have lived than if there will be 
  many trillion persons. Therefore, by Bayes' theorem, you should 
  update your beliefs about mankind's prospects and realise that an 
  impending doomsday is much more probable than you have hitherto thought."
  """
  
  Cf my Gamble model gamble_doomsday.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.

/*
  Using 10**6 balls in the second urn:

  [ball_obs = 7,num_balls1 = 10,num_balls2 = 1000000]
  var : urn
  Probabilities:
  urn1: 0.9997983464408147
  urn2: 0.0002016535591853
  mean = [urn1 = 0.999798,urn2 = 0.000201654]


  But let's make this a little easier and just place 1000 balls in the 
  second urn:

  [ball_obs = 7,num_balls1 = 10,num_balls2 = 1000]
  var : urn
  Probabilities:
  urn1: 0.9891881265972086
  urn2: 0.0108118734027914
  mean = [urn1 = 0.989188,urn2 = 0.0108119]


  [ball_obs = 8,num_balls1 = 10,num_balls2 = 1000]
  var : urn
  Probabilities:
  urn1: 0.9900406504065040
  urn2: 0.0099593495934959
  mean = [urn1 = 0.990041,urn2 = 0.00995935]


  [ball_obs = 9,num_balls1 = 10,num_balls2 = 1000]
  var : urn
  Probabilities:
  urn1: 0.9896517412935323
  urn2: 0.0103482587064677
  mean = [urn1 = 0.989652,urn2 = 0.0103483]


  [ball_obs = 10,num_balls1 = 10,num_balls2 = 1000]
  var : urn
  Probabilities:
  urn1: 0.9913928012519562
  urn2: 0.0086071987480438
  mean = [urn1 = 0.991393,urn2 = 0.0086072]


  [ball_obs = 11,num_balls1 = 10,num_balls2 = 1000]
  var : urn
  Probabilities:
  urn2: 1.0000000000000000
  mean = [urn2 = 1.0]


*/
go ?=>
  member(BallObs,7..11),
  NumBalls1 = 10,
  NumBalls2 = 10**3, % 10**6,
  println([ball_obs=BallObs,num_balls1=NumBalls1,num_balls2=NumBalls2]),
  reset_store,
  run_model(100_000,$model(BallObs,NumBalls1,NumBalls2),[show_probs_trunc,mean]),
  nl,
  % show_store_lengths,nl,
  fail,
  nl.
go => true.
  
  
model(BallObs,NumBalls1,NumBalls2) =>
  Urn = categorical([1/2,1/2],[urn1,urn2]),
  Ball = cond(Urn == urn1,
              random_integer1(NumBalls1),
              random_integer1(NumBalls2)),
  observe(Ball == BallObs),  
  if observed_ok then
    add("urn",Urn),
  end.