/* 

  Firings in Picat.

  From "Resampling Stats Illustrations" (http://www.statistics101.net/PeterBruce_05-illus.pdf)
  Page 54
  """
  Unusual Statistics
  (program “firings”)
  One of the advantages of resampling is its suitability for use with
  non-standard statistics. Here is an illustration of a statistic de-
  signed to meet the needs of a specific situation:
  A company has been accused of firing workers (it has 50) when
  they get close to the level of seniority at which their pension would
  be vested (25 years). The union notes that the levels of seniority of
  7 fired workers in the last 12 months were unusually close to 25
  years.
  Seniority at discharge (years): 23 19 24 23 25 2 5
  Seniority of all workers:
     11 8 24 36 20 19 11 9 10 9 5
     4 2 1 9 21 16 17 11 1 1 23
     19 24 40 28 5 7 1 34 20 16 31
     23 50 4 1 8 8 14 12 32 1 15
     12 25 19 5 24 2
  Note: A “25” indicates the worker’s pension has vested.
  The company counters that operational considerations were the
  only factors in each of the firings and that the proximity of the
  firing dates to pension vesting dates was purely coincidental, the
  result of random chance.
  Can we assess whether this claim is reasonable?
  ...
  Result:
  prob = .11
  The estimated p-value is .11, indicating that a sum as low as the
  observed value of 79 might happen 11% of the time, simply
  drawing workers at random. We conclude the evidence is not
  strong that there was systematic firing of those close to vesting.
  """

  Here we get a little lower probability, but it's still fairly significant: 7.5%

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

/*

  sum_discharged = 79
  var : s
  Probabilities (truncated):
  112: 0.0204000000000000
  108: 0.0204000000000000
  109: 0.0202000000000000
  107: 0.0196000000000000
  .........
  46: 0.0001000000000000
  45: 0.0001000000000000
  43: 0.0001000000000000
  41: 0.0001000000000000
  mean = 108.921

  var : p
  Probabilities:
  false: 0.9249000000000001
  true: 0.0751000000000000
  mean = [false = 0.9249,true = 0.0751]


*/
go ?=>
  SeniorityDischarged = [23,19,24,23,25,2,5],
  SeniorityAllWorkers = [11,8,24,36,20,19,11,9,10,9,5,
                         4,2,1,9,21,16,17,11,1,1,23,
                         19,24,40,28,5,7,1,34,20,16,31,
                         23,50,4,1,8,8,14,12,32,1,15,
                         12,25,19,5,24,2],

  SumDischarged = sum_seniority(SeniorityDischarged),
  println(sum_discharged=SumDischarged),
  reset_store,
  run_model(10_000,$model(SeniorityDischarged,SeniorityAllWorkers,SumDischarged),
                   [show_probs_trunc,mean]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.


sum_seniority(Data) = [cond(25 - V <= 0, 25, 25-V)   : V in Data].sum.
  
model(SeniorityDischarged,SeniorityAllWorkers,SumDischarged) =>
  DischargedSample = draw_without_replacement(SeniorityDischarged.len,
                                              SeniorityAllWorkers),
  S = sum_seniority(DischargedSample),
  P = check(S <= SumDischarged),

  add("s",S),
  add("p",P).