/* 

  7 deaths in one month in Picat.

  From Norman Fenton
  https://x.com/profnfenton/status/1963341658038423939
  """
  With reports of multiple AfD candidates (now up to 7) in Germany dying in the 
  weeks before the election, lots of people have asked me what's the probability 
  that could happen purely by chance. 

  Assuming there are about 3,000 AfD candidates - mostly 'middle-aged' and that 
  7 died within 4-weeks, we have to calculate the probability that at least 7 out of 
  3,000 randomly selected middle-aged Germans would die within a specific 4-week period. 
  This involves some tricky Binomial calculations, so as I'm pressed for time I decided to 
  put the question to ChatGPT (I've never aksed it a probability question before) and was 
  surprised how good a job it did (it used the same method I would have done).

  ... 

  [From ChatGPT:
   If we take a typocal midlife all-cause annual mortality of ~0.5% per year (0.005)
   and scale it to 4 weeks (=28/365) then
    p = 0.005 * 28/365 ~ 0.0000384
  ]
  """
  


  Cf my Gamble model gamble_7_deaths_in_one_month.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, ppl_common_utils.
import util.
% import ordset.

main => go.

/*
  annual_mortality = 0.005
  garbage_collect_when = 100
  var : total deaths
  Probabilities:
  1: 0.3633000000000000
  0: 0.3183000000000000
  2: 0.2080800000000000
  3: 0.0806100000000000
  4: 0.0228900000000000
  5: 0.0056900000000000
  6: 0.0009700000000000
  7: 0.0001400000000000
  9: 0.0000100000000000
  8: 0.0000100000000000
  mean = 1.14827
  HPD intervals:
  HPD interval (0.84): 0.00000000000000..2.00000000000000

  var : p
  Probabilities:
  false: 0.9999800000000000
  true: 0.0000200000000000
  mean = 0.00002
  HPD intervals:
  []

  annual_mortality = 0.007
  garbage_collect_when = 100
  var : total deaths
  Probabilities:
  1: 0.3211800000000000
  2: 0.2616300000000000
  0: 0.1977100000000000
  3: 0.1384100000000000
  4: 0.0568500000000000
  5: 0.0178500000000000
  6: 0.0048400000000000
  7: 0.0012600000000000
  8: 0.0002400000000000
  9: 0.0000300000000000
  mean = 1.61637
  HPD intervals:
  HPD interval (0.84): 0.00000000000000..3.00000000000000

  var : p
  Probabilities:
  false: 0.9997300000000000
  true: 0.0002700000000000
  mean = 0.00027
  HPD intervals:
  []

  annual_mortality = 0.01
  garbage_collect_when = 100
  var : total deaths
  Probabilities (truncated):
  2: 0.2675000000000000
  1: 0.2306500000000000
  3: 0.2031700000000000
  4: 0.1163400000000000
  .........
  8: 0.0016500000000000
  9: 0.0005000000000000
  10: 0.0001200000000000
  11: 0.0000300000000000
  mean = 2.29264
  HPD intervals:
  HPD interval (0.84): 0.00000000000000..4.00000000000000

  var : p
  Probabilities:
  false: 0.9977000000000000
  true: 0.0023000000000000
  mean = 0.0023
  HPD intervals:
  []



*/
go ?=>
  AnnualMortality = 0.005,
  % AnnualMortality = 0.007,
  % AnnualMortality = 0.01,
  println(annual_mortality=AnnualMortality),
  reset_store,
  run_model(100_000,$model(AnnualMortality),[show_probs_trunc,mean,
                           show_hpd_intervals,hpd_intervals=[0.84]
                          ]),
  nl,
  nl.
go => true.
  
  
model(AnnualMortality) =>
  TotalPeople = 3_000,
  N = 7,

  Mortality = 28/365 * AnnualMortality,

  People = bern_n(Mortality,TotalPeople),
  TotalDeaths = People.sum,
  P = check(TotalDeaths > N),
  
  add("total deaths",TotalDeaths),
  add("p",P).