/* 

  A/B testing, given lists of successes for A and B in Picat.

  Cf my Gamble model gamble_ab_testing2.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 : pa
  Probabilities (truncated):
  0.16008925794176: 0.0100000000000000
  0.136774256562123: 0.0100000000000000
  0.133205854485125: 0.0100000000000000
  0.130969519406847: 0.0100000000000000
  .........
  0.021954165962305: 0.0100000000000000
  0.020507910940593: 0.0100000000000000
  0.016715744946103: 0.0100000000000000
  0.008094155321053: 0.0100000000000000
  mean = 0.075073
  [len = 100,min = 0.00809416,mean = 0.075073,median = 0.0768238,max = 0.160089,variance = 0.00127496,stdev = 0.0357066]
  HPD intervals:
  HPD interval (0.94): 0.02195416596230..0.13096951940685

  var : pb
  Probabilities (truncated):
  0.09325977081372: 0.0100000000000000
  0.092415906298853: 0.0100000000000000
  0.088476596564178: 0.0100000000000000
  0.086979360407659: 0.0100000000000000
  .........
  0.00944497182133: 0.0100000000000000
  0.009104191169496: 0.0100000000000000
  0.008684739809529: 0.0100000000000000
  0.00794305794147: 0.0100000000000000
  mean = 0.0429378
  [len = 100,min = 0.00794306,mean = 0.0429378,median = 0.0432625,max = 0.0932598,variance = 0.000467148,stdev = 0.0216136]
  HPD intervals:
  HPD interval (0.94): 0.00794305794147..0.07912587938066

  var : diff
  Probabilities (truncated):
  0.114209696543155: 0.0100000000000000
  0.113817227787623: 0.0100000000000000
  0.100020608360596: 0.0100000000000000
  0.093425781198686: 0.0100000000000000
  .........
  -0.044427200544738: 0.0100000000000000
  -0.056572645717971: 0.0100000000000000
  -0.060275951727162: 0.0100000000000000
  -0.067944386767587: 0.0100000000000000
  mean = 0.0321352
  [len = 100,min = -0.0679444,mean = 0.0321352,median = 0.0296942,max = 0.11421,variance = 0.0018993,stdev = 0.043581]
  HPD intervals:
  HPD interval (0.94): -0.04442720054474..0.09342578119869

  var : a > b
  Probabilities:
  true: 0.7600000000000000
  false: 0.2400000000000000
  mean = 0.76
  show_simple_stats: no numeric data
  HPD intervals:
  show_hpd_intervals: data is not numeric

*/
go ?=>
  %  Generated by ./gener_data.pl
  % True prob 0.05 
  % random_binomial(750,1,0.05);   
  As = [0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,1,0,0,
0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0],
  
  %  Generated by ./gener_data.pl
  % True prob = 0.04
  % random_binomial(750,1,0.04);
  Bs = [0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,1,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0],

  reset_store,
  run_model(10_000,$model(As,Bs),[show_probs_trunc,mean,
                           show_simple_stats,
                           % show_percentiles,
                           show_hpd_intervals,hpd_intervals=[0.94],
                           % show_histogram,
                           min_accepted_samples=100,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model(As,Bs) =>
  PA = beta_dist(1,1), % true prob 0.05
  PB = beta_dist(1,1), % true prob 0.04

  A = bern_n(PA,As.len),
  observe_abc(As,A,1/2),

  B = bern_n(PB,Bs.len),
  observe_abc(Bs,B,1/2),

  Diff = PA - PB,
  AgtB = check(Diff > 0),
  
  if observed_ok then
    add("pa",PA),
    add("pb",PB),
    add("diff",Diff),
    add("a > b",AgtB),    
  end.