/* 

  Battery comparison (Resampling Stats) in Picat.

  From "Resampling Stats Illustrations" (http://www.statistics101.net/PeterBruce_05-illus.pdf)
  Page 40
  """
  Hypothesis test for a difference in means
  (program “battery”)
  Does one brand of car battery last longer than another? Here are
  the figures on 10 randomly selected batteries each for brand A and
  brand B. Is brand A’s apparent advantage significant?
  TABLE 2. BATTERY LIFE TIMES
  brand life (months) aveage
  A 30 32 31 28 31 29 29 24 30 31 29.5
  B 28 28 32 31 24 23 31 27 27 31 28.2
  A’s advantage: 1.3 months
  ...
  We calculate the average life time of each group, and determine whether 
  they differ by as much as the observed data.
  ->
  prob = .147
  """

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

main => go.

/*

  var : a sample mean
  Probabilities (truncated):
  29.681818181818183: 0.0238000000000000
  29.65909090909091: 0.0231000000000000
  29.613636363636363: 0.0221000000000000
  29.477272727272727: 0.0216000000000000
  .........
  27.931818181818183: 0.0001000000000000
  27.84090909090909: 0.0001000000000000
  27.795454545454547: 0.0001000000000000
  27.431818181818183: 0.0001000000000000
  mean = 29.5018
  HPD intervals:
  HPD interval (0.95): 28.63636363636364..30.31818181818182

  var : b sample mean
  Probabilities (truncated):
  28.418181818181822: 0.0091000000000000
  28.145454545454548: 0.0089000000000000
  28.190909090909091: 0.0088000000000000
  27.872727272727275: 0.0086000000000000
  .........
  26.100000000000001: 0.0001000000000000
  26.045454545454547: 0.0001000000000000
  26.027272727272724: 0.0001000000000000
  25.763636363636362: 0.0001000000000000
  mean = 28.1926
  HPD intervals:
  HPD interval (0.95): 27.00909090909091..29.33636363636364

  var : diff
  Probabilities (truncated):
  1.263636363636362: 0.0043000000000000
  1.354545454545452: 0.0033000000000000
  1.127272727272725: 0.0033000000000000
  1.536363636363635: 0.0032000000000000
  .........
  -1.145454545454548: 0.0001000000000000
  -1.154545454545456: 0.0001000000000000
  -1.254545454545454: 0.0001000000000000
  -1.663636363636364: 0.0001000000000000
  mean = 1.30917
  HPD intervals:
  HPD interval (0.95): -0.13636363636364..2.71818181818182

  var : p
  Probabilities:
  1: 0.5071000000000000
  0: 0.4929000000000000
  mean = 0.5071
  HPD intervals:
  HPD interval (0.95): 0.00000000000000..1.00000000000000

  var : p1
  Probabilities:
  1: 0.9622000000000001
  0: 0.0378000000000000
  mean = 0.9622
  HPD intervals:
  HPD interval (0.95): 1.00000000000000..1.00000000000000

  var : p2
  Probabilities:
  0: 0.9631000000000000
  1: 0.0369000000000000
  mean = 0.0369
  HPD intervals:
  HPD interval (0.95): 0.00000000000000..0.00000000000000


*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean,show_hpd_intervals,hpd_intervals=[0.95]]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  TypeA = [30,32,31,28,31,29,29,24,30,31,29.5],
  TypeB = [28,28,32,31,24,23,31,27,27,31,28.2],
  OrigDiff = TypeA.mean - TypeB.mean,
  BothTypes = TypeA ++ TypeB,
  BothTypesLen = BothTypes.len,

  % Draw 10 samples from the combined data
  ASampleMean = resample(BothTypesLen,TypeA).mean,
  BSampleMean = resample(BothTypesLen,TypeB).mean,
  Diff = ASampleMean - BSampleMean,

  P  = check1(Diff >= OrigDiff),
  P1 = check1(ASampleMean > BSampleMean),
  P2 = check1(BSampleMean > ASampleMean),
  
  add("a sample mean",ASampleMean),
  add("b sample mean",BSampleMean),
  add("diff",Diff),
  add("p",P),
  add("p1",P1),
  add("p2",P2).


/*
  Alternative model: use draw_without_replacement instead of resample (see page 42f)

  var : a sample mean
  Probabilities (truncated):
  28.745454545454546: 0.0179000000000000
  28.772727272727273: 0.0174000000000000
  28.836363636363636: 0.0172000000000000
  29.045454545454547: 0.0170000000000000
  .........
  27.272727272727273: 0.0001000000000000
  27.227272727272727: 0.0001000000000000
  27.199999999999999: 0.0001000000000000
  27.109090909090909: 0.0001000000000000
  mean = 28.8462
  HPD intervals:
  HPD interval (0.95): 27.74545454545455..29.88181818181818

  var : b sample mean
  Probabilities (truncated):
  28.954545454545453: 0.0179000000000000
  28.927272727272726: 0.0174000000000000
  28.863636363636363: 0.0172000000000000
  28.654545454545453: 0.0170000000000000
  .........
  27.227272727272727: 0.0001000000000000
  27.181818181818183: 0.0001000000000000
  27.018181818181816: 0.0001000000000000
  26.881818181818179: 0.0001000000000000
  mean = 28.8538
  HPD intervals:
  HPD interval (0.95): 27.77272727272727..29.90909090909091

  var : diff
  Probabilities (truncated):
  -0.209090909090907: 0.0179000000000000
  -0.154545454545453: 0.0174000000000000
  -0.027272727272727: 0.0172000000000000
  0.390909090909094: 0.0170000000000000
  .........
  -3.154545454545453: 0.0001000000000000
  -3.245454545454546: 0.0001000000000000
  -3.300000000000001: 0.0001000000000000
  -3.481818181818181: 0.0001000000000000
  mean = -0.00764545
  HPD intervals:
  HPD interval (0.95): -2.20909090909091..2.06363636363636

  var : p
  Probabilities:
  0: 0.8722000000000000
  1: 0.1278000000000000
  mean = 0.1278
  HPD intervals:
  HPD interval (0.95): 0.00000000000000..1.00000000000000

  var : p1
  Probabilities:
  0: 0.5028000000000000
  1: 0.4972000000000000
  mean = 0.4972
  HPD intervals:
  HPD interval (0.95): 0.00000000000000..1.00000000000000

  var : p2
  Probabilities:
  1: 0.5028000000000000
  0: 0.4972000000000000
  mean = 0.5028
  HPD intervals:
  HPD interval (0.95): 0.00000000000000..1.00000000000000

*/
go2 ?=>
  reset_store,
  run_model(10_000,$model2,[show_probs_trunc,mean,show_hpd_intervals,hpd_intervals=[0.95]]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go2 => true.
  
  
model2() =>
  TypeA = [30,32,31,28,31,29,29,24,30,31,29.5],
  TypeB = [28,28,32,31,24,23,31,27,27,31,28.2],
  OrigDiff = TypeA.mean - TypeB.mean,
  
  BothTypes = TypeA ++ TypeB,

  BothSamples = draw_without_replacement(BothTypes),

  % Pick the first TypeA.len elements
  ASampleMean = take(BothSamples,TypeA.len).mean,
  BSampleMean = drop(BothSamples,TypeA.len).mean,
  Diff = ASampleMean - BSampleMean,

  P  = check1(Diff >= OrigDiff),
  P1 = check1(ASampleMean > BSampleMean),
  P2 = check1(BSampleMean > ASampleMean),
  
  add("a sample mean",ASampleMean),
  add("b sample mean",BSampleMean),
  add("diff",Diff),
  add("p",P),
  add("p1",P1),
  add("p2",P2).