/* 

  Holmes Clock problem in Picat.

  From "Introduction to Probability" av Grinstead & Snell, page 92
  """
  Mr. Wimply Dimple, one of London's most prestigious watch makers, has come to 
  Sherlock Holmes in a panic, having discovered that someone has been producing and 
  selling crude counterfeits of his best selling watch. The 16 counterfeits so far 
  discovered bear stamped numbers, all of which fall between 1 and 56, and Dimple is 
  anxious to know the extent of the forger's work. All present agree that it seems 
  reasonable to assume that the counterfeits thus far produced bear consecutive 
  numbers from 1 to whatever the total number is. "Chin up, Dimple," opines Dr. Watson.
  "I shouldn't worry overly much if I were you; the Maximum Likelihood Principle, 
  which estimates the total number as precisely that which gives the highest 
  probability for the series of numbers found, suggests that we guess 56 itself as the
  total. Thus, your forgers are not a big operation, and we shall have them safely 
  behind bars before your business suffers significantly." "Stuff, nonsense, and 
  bother your fancy principles, Watson," counters Holmes. "Anyone can see that, of 
  course, there must be quite a few more than 56 watches - why the odds of our having 
  discovered precisely the highest numbered watch made are laughably negligible. A 
  much better guess would be twice 56." 
  """

  Let's simulate this by an urn model were we have two urns:
  - urn 1 which have 56 balls numbered 1..56
  - urn 2 which have 2*56 = 112 balls numbered 1..112

  The probability of selecting each of the urns are 1/2.

  Using Bayes formula:

      p(ball from urn 1) = 0.5
      p(ball from urn 2) = 0.5

      p(number 56 is drawn|urn=urn 1) = 1/56
      p(nummer 56 is drawn|urn=urn 2) = 1/112

  What are the probabilities:
    p(urn 1|number 56 was drawn)
    p(urn 2|number 56 was drawn)

  p(urn 1|number 56 was drawn) = p(urn 1 /\ number 56) / p(number 56)

                          = p(number 56|urn 1) * p(urn 1)
                            ------------------------------
                            p(number 56|urn 1) * p(urn 1) + p(number 56|urn 2)*p(urn 2)

                          = 1/56 * 0.5
                            ----------
                            1/56 * 0.5 + 1/112*0.5

                          = (1/56 * 0.5)/(1/56 * 0.5 + 1/112*0.5) 
                          = 0.6666667

  p(urn 2|number 56 was drawn) = p(urn 2 /\ number 56) / p(number 56)

                          = p(number 56|urn 2) * p(urn 2)
                            ------------------------------
                            p(number 56|urn 1) * p(urn 1) + p(number 56|urn 2)*p(urn 2)

                          = 1/112 * 0.5
                            ----------
                            1/56 * 0.5 + 1/112*0.5

                          = (1/112 * 0.5)/(1/56 * 0.5 + 1/112*0.5) 
                          = 0.3333333

  This model gives the following result which confirms the calculations 
  (or the other way around :-):

   var : urn
   Probabilities:
   urn1: 0.6685878962536023
   urn2: 0.3314121037463977
   mean = [urn1 = 0.668588,urn2 = 0.331412]

  I.e. if we draw ball number 56 it's more likely that it was from urn number 1 
  (with balls 1..56) than from urn 2 (with balls 1..112), 66% vs 33%.


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

/*

  Here are some more experiments were we increments the number of balls 
  in urn2:

  [numballs1 = 56,ball_obs = 56]
  num_balls2 = 112
  var : urn
  Probabilities:
  urn1: 0.6543396226415095
  urn2: 0.3456603773584906
  mean = [urn1 = 0.65434,urn2 = 0.34566]

  num_balls2 = 168
  var : urn
  Probabilities:
  urn1: 0.7433198380566801
  urn2: 0.2566801619433198
  mean = [urn1 = 0.74332,urn2 = 0.25668]

  num_balls2 = 224
  var : urn
  Probabilities:
  urn1: 0.8018181818181818
  urn2: 0.1981818181818182
  mean = [urn1 = 0.801818,urn2 = 0.198182]

  num_balls2 = 280
  var : urn
  Probabilities:
  urn1: 0.8518867924528302
  urn2: 0.1481132075471698
  mean = [urn1 = 0.851887,urn2 = 0.148113]

  num_balls2 = 336
  var : urn
  Probabilities:
  urn1: 0.8496605237633366
  urn2: 0.1503394762366634
  mean = [urn1 = 0.849661,urn2 = 0.150339]

  num_balls2 = 392
  var : urn
  Probabilities:
  urn1: 0.8821782178217822
  urn2: 0.1178217821782178
  mean = [urn1 = 0.882178,urn2 = 0.117822]

  num_balls2 = 448
  var : urn
  Probabilities:
  urn1: 0.8841463414634146
  urn2: 0.1158536585365854
  mean = [urn1 = 0.884146,urn2 = 0.115854]

  num_balls2 = 504
  var : urn
  Probabilities:
  urn1: 0.9056203605514316
  urn2: 0.0943796394485684
  mean = [urn1 = 0.90562,urn2 = 0.0943796]

  num_balls2 = 560
  var : urn
  Probabilities:
  urn1: 0.9229188078108942
  urn2: 0.0770811921891059
  mean = [urn1 = 0.922919,urn2 = 0.0770812]

*/
go ?=>
  NumBalls1 = 56,
  BallObs = 56,
  println([numballs1=NumBalls1,ball_obs=56]),
  member(NumBalls2,[112,168,224,280,336,392,448,504,560]),
  println(num_balls2=NumBalls2),
  reset_store,
  run_model(100_000,$model(NumBalls1,NumBalls2,BallObs),[show_probs_trunc,mean]),
  nl,
  % show_store_lengths,
  fail,
  nl.
go => true.
  
  
model(NumBalls1,NumBalls2,BallObs) =>
  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.