/* 

  IQ best in Picat.

  Probability that the best male IQ is better than the best female IQ,
  given that 
  - male IQ is normal 100 15 
  - female IQ is normal 100 10

  Cf 
  - ppl_spread_of_traits.pi
  - my Gamble model gamble_iq_best.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.

/*
  For some different sizes of population

  n = 10
  var : best man
  Probabilities (truncated):
  163.496224469514004: 0.0001000000000000
  161.543889634799768: 0.0001000000000000
  159.996516898252764: 0.0001000000000000
  159.447481673550328: 0.0001000000000000
  .........
  97.532609977468525: 0.0001000000000000
  96.886377189744948: 0.0001000000000000
  96.882959162449424: 0.0001000000000000
  95.528336257117232: 0.0001000000000000
  mean = 123.124
  [len = 10000,min = 95.5283,mean = 123.124,median = 122.587,max = 163.496,variance = 76.8187,stdev = 8.76463]
  HPD intervals:
  HPD interval (0.84): 110.12186680317701..134.44796680952743
  HPD interval (0.9): 108.26184769265585..136.72975582859686
  HPD interval (0.94): 107.45978075807147..139.87833876751316
  HPD interval (0.99): 103.52974056637925..148.11545234111156
  HPD interval (0.99999): 95.52833625711723..163.49622446951400

  var : best woman
  Probabilities (truncated):
  144.351552254739374: 0.0001000000000000
  143.326803932325674: 0.0001000000000000
  140.268501740491416: 0.0001000000000000
  138.949514766520849: 0.0001000000000000
  .........
  100.344072374181309: 0.0001000000000000
  100.129730212294902: 0.0001000000000000
  99.683600407033765: 0.0001000000000000
  98.456465976706951: 0.0001000000000000
  mean = 115.401
  [len = 10000,min = 98.4565,mean = 115.401,median = 114.965,max = 144.352,variance = 33.9501,stdev = 5.82667]
  HPD intervals:
  HPD interval (0.84): 106.96962333127327..123.12338501733680
  HPD interval (0.9): 105.59995782264727..124.54593005938723
  HPD interval (0.94): 105.17777091493218..126.74157817286873
  HPD interval (0.99): 101.88822397345699..131.46608461478820
  HPD interval (0.99999): 98.45646597670695..144.35155225473937

  var : p1
  Probabilities:
  true: 0.7660000000000000
  false: 0.2340000000000000
  mean = 0.766

  var : p2
  Probabilities:
  false: 0.7660000000000000
  true: 0.2340000000000000
  mean = 0.234


  n = 100
  var : best man
  Probabilities (truncated):
  170.160992361526723: 0.0001000000000000
  168.267420966701565: 0.0001000000000000
  166.512677917074541: 0.0001000000000000
  166.246446729366994: 0.0001000000000000
  .........
  121.3931833098152: 0.0001000000000000
  121.357570213388243: 0.0001000000000000
  120.954359813686224: 0.0001000000000000
  119.566316565726567: 0.0001000000000000
  mean = 137.57
  [len = 10000,min = 119.566,mean = 137.57,median = 136.852,max = 170.161,variance = 40.9882,stdev = 6.4022]
  HPD intervals:
  HPD interval (0.84): 128.28496728550996..145.49066495352071
  HPD interval (0.9): 127.65369236143340..147.93931824580494
  HPD interval (0.94): 125.85586990394341..149.20362693164697
  HPD interval (0.99): 124.07425650121847..156.63815830481178
  HPD interval (0.99999): 119.56631656572657..170.16099236152672

  var : best woman
  Probabilities (truncated):
  149.87028677536847: 0.0001000000000000
  146.197696144211193: 0.0001000000000000
  143.584072027420746: 0.0001000000000000
  143.410368817326685: 0.0001000000000000
  .........
  113.604381465534473: 0.0001000000000000
  113.554415233791673: 0.0001000000000000
  113.427924594947427: 0.0001000000000000
  112.838627196865261: 0.0001000000000000
  mean = 124.977
  [len = 10000,min = 112.839,mean = 124.977,median = 124.506,max = 149.87,variance = 18.8991,stdev = 4.34731]
  HPD intervals:
  HPD interval (0.84): 118.89058287473236..130.48843394812641
  HPD interval (0.9): 118.35384646109104..132.18988574056428
  HPD interval (0.94): 117.08134233114021..133.12872002515908
  HPD interval (0.99): 115.49998134831770..138.04884636901519
  HPD interval (0.99999): 112.83862719686526..149.87028677536847

  var : p1
  Probabilities:
  true: 0.9567000000000000
  false: 0.0433000000000000
  mean = 0.9567

  var : p2
  Probabilities:
  false: 0.9567000000000000
  true: 0.0433000000000000
  mean = 0.0433


  n = 1000
  var : best man
  Probabilities (truncated):
  178.566660324808538: 0.0001000000000000
  175.368022573924975: 0.0001000000000000
  173.606947129934667: 0.0001000000000000
  172.244321130071398: 0.0001000000000000
  .........
  137.059568757735377: 0.0001000000000000
  136.994220626600907: 0.0001000000000000
  136.936961693364395: 0.0001000000000000
  135.682306840853954: 0.0001000000000000
  mean = 148.53
  [len = 10000,min = 135.682,mean = 148.53,median = 147.954,max = 178.567,variance = 26.3073,stdev = 5.12906]
  HPD intervals:
  HPD interval (0.84): 141.22525370165368..155.03639539203888
  HPD interval (0.9): 140.36792118873359..156.63032601229841
  HPD interval (0.94): 139.51514035243102..158.17326125534211
  HPD interval (0.99): 137.72927979782162..163.78177402953756
  HPD interval (0.99999): 135.68230684085395..178.56666032480854

  var : best woman
  Probabilities (truncated):
  151.09047321790564: 0.0001000000000000
  150.586072802181064: 0.0001000000000000
  149.286722679755513: 0.0001000000000000
  149.213979617441936: 0.0001000000000000
  .........
  124.291600401669314: 0.0001000000000000
  124.174379629763848: 0.0001000000000000
  124.101200843669076: 0.0001000000000000
  123.7651945110171: 0.0001000000000000
  mean = 132.398
  [len = 10000,min = 123.765,mean = 132.398,median = 131.923,max = 151.09,variance = 12.6624,stdev = 3.55843]
  HPD intervals:
  HPD interval (0.84): 127.49126728740380..136.77651041536444
  HPD interval (0.9): 126.86093234722195..137.94618603243686
  HPD interval (0.94): 126.43799078495221..139.20959874032712
  HPD interval (0.99): 125.51774835189796..144.03975889818526
  HPD interval (0.99999): 123.76519451101710..151.09047321790564

  var : p1
  Probabilities:
  true: 0.9957000000000000
  false: 0.0043000000000000
  mean = 0.9957

  var : p2
  Probabilities:
  false: 0.9957000000000000
  true: 0.0043000000000000
  mean = 0.0043


*/
go ?=>
  member(N, [10,100,1000]),
  println(n=N),
  reset_store,
  run_model(10_000,$model(N),[show_probs_trunc,mean,
                           show_simple_stats,
                           % show_percentiles,
                           show_hpd_intervals,hpd_intervals=[0.84,0.9,0.94,0.99,0.99999]
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  fail,
  nl.
go => true.
  
  
model(N) =>
  Men = normal_dist_n(100,15,N),
  Women = normal_dist_n(100,10,N),

  BestMan = Men.max,
  BestWoman = Women.max,
  P1 = check(BestMan > BestWoman),
  P2 = check(BestWoman > BestMan),  
  
  add("best man",BestMan),
  add("best woman",BestWoman),
  add("p1",P1),
  add("p2",P2).