/* 

  All Girls world puzzle in Picat.

  From https://brainstellar.com/puzzles/probability/6
  """
  All Girls World?

  In a world where everyone wants a girl child, each family continues having 
  babies till they have a girl. What do you think will the boy-to-girl ratio 
  be eventually?

  Assuming probability of having a boy or a girl is the same and there is no 
  other gender at the time of birth.

  
  Answer:
  Suppose there are N couples. First time, N/2 girls and N/2 boys are born. N/2 
  couples retire, and rest half try for another child.

  Next time, N/4 couples give birth to N/4 girls and N/4 boys. Thus, even in the 
  second iteration, the ratio is 1:1. It can now be seen that this ratio will 
  always remain the same, no matter how many times people try to give birth to a 
  favored gender.
  """

  This is a port of my Gamble model gamble_all_girls_world.rkt which also includes
  more exact probabilities.

  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.

/*

  var : len
  Probabilities:
  1: 0.50053333333333339 (1877 / 3750)
  2: 0.24990000000000001 (2499 / 10000)
  3: 0.12220000000000000 (611 / 5000)
  4: 0.06130000000000000 (613 / 10000)
  5: 0.03310000000000000 (331 / 10000)
  6: 0.01676666666666667 (503 / 30000)
  7: 0.00880000000000000 (11 / 1250)
  8: 0.00376666666666667 (113 / 30000)
  9: 0.00180000000000000 (9 / 5000)
  10: 0.00096666666666667 (29 / 30000)
  11: 0.00043333333333333 (13 / 30000)
  12: 0.00020000000000000 (1 / 5000)
  14: 0.00010000000000000 (1 / 10000)
  15: 0.00006666666666667 (1 / 15000)
  13: 0.00006666666666667 (1 / 15000)
  mean = 2.00627

  var : ratio
  Probabilities:
  0: 0.50053333333333339 (1877 / 3750)
  1: 0.24990000000000001 (2499 / 10000)
  2: 0.12220000000000000 (611 / 5000)
  3: 0.06130000000000000 (613 / 10000)
  4: 0.03310000000000000 (331 / 10000)
  5: 0.01676666666666667 (503 / 30000)
  6: 0.00880000000000000 (11 / 1250)
  7: 0.00376666666666667 (113 / 30000)
  8: 0.00180000000000000 (9 / 5000)
  9: 0.00096666666666667 (29 / 30000)
  10: 0.00043333333333333 (13 / 30000)
  11: 0.00020000000000000 (1 / 5000)
  13: 0.00010000000000000 (1 / 10000)
  14: 0.00006666666666667 (1 / 15000)
  12: 0.00006666666666667 (1 / 15000)
  mean = 1.00627


  The ratio is thus ~ 1 (with a possible limit of 20 children per family).

  Cf the geometric distribution:
  Picat> X=[1-geometric_dist_cdf(1/2,I) :  I in 0..10]
X = [0.5,0.25,0.125,0.0625,0.03125,0.015625,0.0078125,0.00390625,0.001953125,0.0009765625,0.00048828125]

  A more exact probability:

  Picat> X=[1-geometric_dist_cdf(1/2,I) :  I in 0..10].sum
  X = 0.99951171875

  Picat> X=[1-geometric_dist_cdf(1/2,I) :  I in 0..100].sum
  X = 1.0

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_rat,mean]),
  % fail,
  nl.
go => true.

family(L) = Ret =>
  C = uniform_draw(["boy","girl"]),
  L2 = L ++ [C],
  if C == "girl" then
    Ret = L2
  else
    Ret = family(L2)
  end.
    
model() =>
  ThisFamily = family([]),
  Len = ThisFamily.len,
  Ratio = Len - 1, % i.e. there are Len - 1 boys and 1 girl
  add("len",Len),
  add("ratio",Ratio).


/*
  Let's run the model a coupl of time to just see the spread of the result.

  ratio = 1.0254
  ratio = 1.0054
  ratio = 1.0019
  ratio = 1.0042
  ratio = 0.9898
  ratio = 0.997
  ratio = 0.9992
  ratio = 0.9846
  ratio = 1.0197
  ratio = 1.0194
  ratio = 0.9884
  ratio = 0.9775
  ratio = 1.0176
  ratio = 1.0286
  ratio = 0.9891
  ratio = 1.004
  ratio = 0.9689
  ratio = 0.9814
  ratio = 0.9847
  ratio = 0.9974
  ratios = [0.9689,0.9775,0.9814,0.9846,0.9847,0.9884,0.9891,0.9898,0.997,0.9974,0.9992,1.0019,1.004,1.0042,1.0054,1.0176,1.0194,1.0197,1.0254,1.0286]
  [len = 20,min = 0.9689,mean = 0.99921,median = 0.9983,max = 1.0286,variance = 0.000262787,stdev = 0.0162107]
  Scatter plot:

  ymax = 1.0286
  |                                                           *
  |                                                        *   
  |                                                            
  |                                                  *  *      
  |                                               *            
  |                                                            
  |                                                            
  |                                           *                
  |                                  *  *  *                   
  |                               *                            
  |                         *  *                               
  |                                                            
  |                      *                                     
  |                *  *                                        
  |         *  *                                               
  |      *                                                     
  |   *                                                        
  |                                                            
  |                                                            
  |*                                                           
  ymin = 0.9689
  x = [1,20]


*/
go2 =>
  Ratios = [],
  foreach(_ in 1..20)
    reset_store,
    run_model(10_000,$model,[mean,presentation=[]]),
    Ratio = get_store().get("ratio").mean,
    println(ratio=Ratio),
    Ratios := Ratios ++ [Ratio],
  end,
  Sorted = Ratios.sort,
  println(ratios=Sorted),
  Sorted.show_simple_stats,
  Sorted.show_scatter_plot,
  nl.