/* 

  Locomotive problem in Picat.

  From Think Bayes, page 22ff
  """
  I found the locomotive problem in Frederick Mosteller's "Fifty Challenging
  Problems in Probability with Solutions" (Dover, 1987):
     'A railroad numbers its locomotives in order 1..N. One day you
      see a locomotive with the number 60. Estimate how many loco-
      motives the railroad has.'

  ...

  The mean of the posterior is 333, so that might be a good guess if you wanted to
  minimize error.
  """
  
  As the book (Think Bayes) mentions, note that it's very sensitive to 
  the value of maxInt. Here are some results for different maxInt 
  (from the book):
  maxInt  N
   -------------------------
    100    77.8166306057923
    200   115.84577808279282
    500   207.2393675826458
   1000   334.04414386751915 
   2000   553.5237331955558

  Mosteller's solution: 2*(60-1)+1: 119.


  This is a port of my Racket/Gamble model gamble_locomotive_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.

main => go.


/*

  mosteller = 119

  maxInt = 100
  Var n:
  Probabilities:
  97: 0.0625000000000000
  90: 0.0625000000000000
  80: 0.0625000000000000
  79: 0.0625000000000000
  77: 0.0625000000000000
  73: 0.0625000000000000
  63: 0.0625000000000000
  61: 0.0625000000000000
  98: 0.0312500000000000
  94: 0.0312500000000000
  91: 0.0312500000000000
  89: 0.0312500000000000
  87: 0.0312500000000000
  83: 0.0312500000000000
  75: 0.0312500000000000
  74: 0.0312500000000000
  72: 0.0312500000000000
  69: 0.0312500000000000
  68: 0.0312500000000000
  66: 0.0312500000000000
  65: 0.0312500000000000
  64: 0.0312500000000000
  62: 0.0312500000000000
  60: 0.0312500000000000
  mean = 76.78125



  maxInt = 200
  Var n:
  Probabilities:
  64: 0.0952380952380952
  196: 0.0476190476190476
  183: 0.0476190476190476
  182: 0.0476190476190476
  179: 0.0476190476190476
  174: 0.0476190476190476
  164: 0.0476190476190476
  149: 0.0476190476190476
  146: 0.0476190476190476
  141: 0.0476190476190476
  132: 0.0476190476190476
  121: 0.0476190476190476
  105: 0.0476190476190476
  104: 0.0476190476190476
  103: 0.0476190476190476
  83: 0.0476190476190476
  78: 0.0476190476190476
  75: 0.0476190476190476
  68: 0.0476190476190476
  60: 0.0476190476190476
  mean = 122.429

  maxInt = 500
  Var n:
  Probabilities:
  485: 0.0476190476190476
  473: 0.0476190476190476
  418: 0.0476190476190476
  410: 0.0476190476190476
  270: 0.0476190476190476
  255: 0.0476190476190476
  220: 0.0476190476190476
  190: 0.0476190476190476
  154: 0.0476190476190476
  143: 0.0476190476190476
  141: 0.0476190476190476
  126: 0.0476190476190476
  120: 0.0476190476190476
  114: 0.0476190476190476
  112: 0.0476190476190476
  111: 0.0476190476190476
  100: 0.0476190476190476
  91: 0.0476190476190476
  74: 0.0476190476190476
  69: 0.0476190476190476
  60: 0.0476190476190476
  mean = 196.952

  maxInt = 1000
  Var n:
  Probabilities:
  981: 0.1111111111111111
  866: 0.1111111111111111
  379: 0.1111111111111111
  259: 0.1111111111111111
  230: 0.1111111111111111
  198: 0.1111111111111111
  193: 0.1111111111111111
  184: 0.1111111111111111
  78: 0.1111111111111111
  mean = 374.222

  maxInt = 2000
  Var n:
  Probabilities:
  1699: 0.1428571428571428
  1511: 0.1428571428571428
  1224: 0.1428571428571428
  1190: 0.1428571428571428
  954: 0.1428571428571428
  632: 0.1428571428571428
  464: 0.1428571428571428
  mean = 1096.29

*/
go ?=>
  Obs = 60,
  println(mosteller=(2*(Obs-1)+1)),
  member(MaxInt,[100,200,500,1000,2000]),
  nl,
  println(maxInt=MaxInt),
  reset_store,
  run_model(30_000,$model(Obs,MaxInt),[show_probs,mean]),
  fail,
  nl.
go => true.

model(Obs,MaxInt) =>
  Priors = [1 / (I+1) : I in 0..MaxInt],
  N = categorical(Priors,0..MaxInt),
  Locomotive = random_integer1(MaxInt),

  observe(Locomotive == Obs, N >= Locomotive),

  if observed_ok  then
    add("n",N)
  end.