/* 

  Comparing two proportions in Picat.

  From http://www.statistics101.net/statistics101web_000007.htm
  """
  From CliffsQuickReview Statistics, Example 16, Page 92: 
  A swimming school wants to determine whether a recently 
  hired instructor is working out. Sixteen out of 25 of 
  Instructor A's students passed the lifeguard certification 
  test on the first try. In comparison, 57 out of 72 of more 
  experienced Instructor B's students passed the test on the 
  first try. Is Instructor A's success rate worse than 
  Instructor B's? Use alpha = 0.10. 
 
  Null hypothesis: A  s rate is >= B  s rate 
  Alternate hypothesis: A  s rate is < B  s rate 
  This is a one-tailed test. 
  """

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

/*

  var : SampleA
  Probabilities (truncated):
  64.0: 0.1629000000000000
  68.0: 0.1587000000000000
  60.0: 0.1481000000000000
  72.0: 0.1235000000000000
  .........
  32.0: 0.0019000000000000
  92.0: 0.0008000000000000
  28.000000000000004: 0.0003000000000000
  96.0: 0.0001000000000000
  mean = 63.9268

  var : SampleB
  Probabilities (truncated):
  80.555555555555557: 0.1193000000000000
  79.166666666666657: 0.1160000000000000
  81.944444444444443: 0.1067000000000000
  77.777777777777786: 0.1060000000000000
  .........
  93.055555555555557: 0.0007000000000000
  62.5: 0.0007000000000000
  94.444444444444443: 0.0001000000000000
  59.722222222222221: 0.0001000000000000
  mean = 79.1469

  var : p
  Probabilities:
  false: 0.9275000000000000
  true: 0.0725000000000000
  mean = [false = 0.9275,true = 0.0725]

*/
go ?=>
  A = ones(16,passed) ++ ones(9,failed),
  B = ones(57,passed) ++ ones(15,failed),    
  reset_store,
  run_model(10_000,$model(A,B),[show_probs_trunc,mean]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model(A,B) =>
  SampleA = 100 * (count_occurrences(passed, resample(25,A)) / 25),
  SampleB = 100 * (count_occurrences(passed, resample(72,B)) / 72),
  P = check(SampleA >= SampleB),
  
  add("SampleA",SampleA),
  add("SampleB",SampleB),
  add("p",P).



/*

  Here we look for the the rates for a and b in binomial(p,n),
  and check the difference of the two rates.

  var : rate a
  Probabilities (truncated):
  0.793970941661941: 0.0217391304347826
  0.779805961238646: 0.0217391304347826
  0.770529913866837: 0.0217391304347826
  0.746316577676847: 0.0217391304347826
  .........
  0.476579435027787: 0.0217391304347826
  0.476366114697649: 0.0217391304347826
  0.46836820214564: 0.0217391304347826
  0.433859351563082: 0.0217391304347826
  mean = 0.610815

  var : rate b
  Probabilities (truncated):
  0.887103468277631: 0.0217391304347826
  0.880571071897608: 0.0217391304347826
  0.880301530013749: 0.0217391304347826
  0.864609184061101: 0.0217391304347826
  .........
  0.733550319183284: 0.0217391304347826
  0.715219662761071: 0.0217391304347826
  0.700195820382628: 0.0217391304347826
  0.695540286585699: 0.0217391304347826
  mean = 0.791902

  var : p
  Probabilities:
  false: 0.9347826086956522
  true: 0.0652173913043478
  mean = [false = 0.934783,true = 0.0652174]

*/
go2 ?=>
  reset_store,
  run_model(100_000,$model2,[show_probs_trunc,mean]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go2 => true.
  
  
model2() =>
  NA = 25, % Number of trials for A
  NB = 72, % Number of trials for A

  RateA = beta_dist(1,1),
  SA = binomial_dist(NA,RateA),
  observe(SA == 16), % number of successes for A

  RateB = beta_dist(1,1),
  SB = binomial_dist(NB,RateB),
  observe(SB == 57), % number of successes for B

  P = check(RateA > RateB),

  if observed_ok then
    add("rate a",RateA),
    add("rate b",RateB),
    add("p",P)
  end.