#! /usr/bin/perl -w
use strict;
use Statistics::Distributions qw(chisqrprob udistr);

my ($n, $trials) = @ARGV;
$n ||= 2;
$trials ||= 10_000;

for (1..$trials) {
  # Simulate an A/B test with n identical versions
  my @measurements
    = map {
        # This is the normal approximation to 1,000,000 random coin flips
        500_000 + 500 * udistr(rand(1));
      } 1..$n;;
  
  # Find the best and the worst.
  my $min = my $max = $measurements[0];
  for my $measurement (@measurements) {
    if ($measurement < $min) {
      $min = $measurement;
    }
    elsif ($measurement > $max) {
      $max = $measurement;
    }
  }

  # Now calculate the g-test for the best and worst
  my $a = my $b = 1_000_000;
  my $a_yes = $max;
  my $a_no = $a - $max;
  my $b_yes = $min;
  my $b_no = $b - $min;

  my $total = $a + $b;
  my $yes = $min + $max;
  my $no = $total - $yes;

  my $e_a_yes = $a * $yes / $total;
  my $e_a_no  = $a * $no  / $total;
  my $e_b_yes = $b * $yes / $total;
  my $e_b_no  = $b * $no  / $total;

  my $g_test =
    2 * (
        $a_yes * log($a_yes / $e_a_yes)
      + $a_no  * log($a_no  / $e_a_no)
      + $b_yes * log($b_yes / $e_b_yes)
      + $b_no  * log($b_no  / $e_b_no)
    );

  my $p = chisqrprob(1, $g_test);

  printf "%.7f\n", $p;
}