{{ message }}

Instantly share code, notes, and snippets.

# stucchio/bayesian_ab_test.py

Last active Nov 29, 2021
Bayesian A/B test code
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 from matplotlib import use from pylab import * from scipy.stats import beta, norm, uniform from random import random from numpy import * import numpy as np import os # Input data prior_params = [ (1, 1), (1,1) ] threshold_of_caring = 0.001 N = array([ 200, 204 ]) s = array([ 16, 36 ]) #Don't edit anything past here def bayesian_expected_error(N,s): degrees_of_freedom = len(prior_params) posteriors = [] for i in range(degrees_of_freedom): posteriors.append( beta(prior_params[i] + s[i] - 1, prior_params[i] + N[i] - s[i] - 1) ) xgrid_size = 1024 x = mgrid[0:xgrid_size,0:xgrid_size] / float(xgrid_size) # Compute joint posterior, which is a product distribution pdf_arr = posteriors.pdf(x) * posteriors.pdf(x) pdf_arr /= pdf_arr.sum() # normalization expected_error_dist = maximum(x-x,0.0) * pdf_arr return expected_error_dist.sum() # Code degrees_of_freedom = len(prior_params) posteriors = [] for i in range(degrees_of_freedom): posteriors.append( beta(prior_params[i] + s[i] - 1, prior_params[i] + N[i] - s[i] - 1) ) if degrees_of_freedom == 2: xgrid_size = 1024 x = mgrid[0:xgrid_size,0:xgrid_size] / float(xgrid_size) # Compute joint posterior, which is a product distribution pdf_arr = posteriors.pdf(x) * posteriors.pdf(x) pdf_arr /= pdf_arr.sum() # normalization prob_error = zeros(shape=x.shape) if (s / float(N)) > (s / float(N)): prob_error[where(x > x)] = 1.0 else: prob_error[where(x > x)] = 1.0 expected_error = maximum(abs(x-x),0.0) expected_err_scalar = (expected_error * prob_error * pdf_arr).sum() if (expected_err_scalar < threshold_of_caring): if (s / float(N)) > (s / float(N)): print "Probability that version B is larger is " + str((prob_error*pdf_arr).sum()) print "Terminate test. Choose version B. Expected error is " + str(expected_err_scalar) else: print "Probability that version A is larger is " + str((prob_error*pdf_arr).sum()) print "Terminate test. Choose version A. Expected error is " + str(expected_err_scalar) else: print "Probability that version B is larger is " + str((prob_error*pdf_arr).sum()) print "Continue test. Expected error was " + str(expected_err_scalar) + " > " + str(threshold_of_caring)

### louispotok commented Aug 15, 2014

Looks like the bayesian_expected_error function is never called. Is there a reason you left it in?

### jjhanlon commented Jan 20, 2015

Large numbers seem to be breaking the code for some reason, e.g. I've entered in the following for N and s

N = array([ 15709030, 1448607 ])
s = array([ 955, 157 ])

This results in:
Probability that version B is larger is nan
Continue test. Expected error was nan > 0.001
/Users/jhanlon/Dropbox/BayesianStat.py:44: RuntimeWarning: invalid value encountered in divide
pdf_arr /= pdf_arr.sum() # normalization

I haven't touched the priors or delta

### arnov commented Feb 13, 2015

I've implemented a python version for this as well. It's a bit simpler, it closely follows Evan Millar's explanation.