Probably buggy A/B test simulation code

ab_simulation.py

#!/usr/bin/env python

import sys
import numpy as np
from scipy.stats import binom, beta, chi2

m     = int(sys.argv[1])      # number of trials
pa    = float(sys.argv[2])    # baseline proportion (b/tw 0 and 1)
t     = float(sys.argv[3])    # materiality threshold (b/tw 0 and 1)
                              # (% increase relative to baseline proportion)
alpha = float(sys.argv[4])    # confidence level
n     = int(sys.argv[5])      # sample size

pb = (1.0 + t) * pa           # increased proportion (goal = detect this)
ea = pa * n                   # expected count for proportion = pa
eb = pa * n                   # expected count for proportion = pb

A  = binom.rvs(n, pa, size=m)
A2 = binom.rvs(n, pa, size=m) # for A/A test
B  = binom.rvs(n, pb, size=m) # for A/B test

x_crit = chi2.isf(alpha, 1)   # critical value, report A =/= B when X^2 > this
X1 = (A - ea)**2 / ea + (A2 - eb)**2 / eb  # A/A test
X2 = (A - ea)**2 / ea +  (B - eb)**2 / eb  # A/B test

freq_I = float(np.sum(X1 > x_crit)) / m * 100.0
freq_II = (1.0 - (float(np.sum(X2 > x_crit)) / m)) * 100.0

# Jeffrey's prior = beta(0.5, 0.5)
# PPF = inverse CDF
A_ub  = beta.ppf(1.0 - alpha, 0.5 + A, 0.5 + n - A) # upper bound for A
A2_lb = beta.ppf(alpha, 0.5 + A2, 0.5 + n - A2)     # lower bound for A2
B_lb  = beta.ppf(alpha, 0.5 + B, 0.5 + n - B)       # lower bound for B

bayes_I = (float(np.sum(A_ub < A2_lb)) / m) * 100.0
bayes_II = (1.0 - (float(np.sum(A_ub < B_lb)) / m)) * 100.0

print "Frequentist  Type I Error    %.1f%%" % freq_I
print "Bayesian     Type I Error    %.1f%%" % bayes_I
print "Frequentist  Type II Error   %.1f%%" % freq_II
print "Bayesian     Type II Error   %.1f%%" % bayes_II