# Bayesian A/B-test having trouble converging in PyMC

I'm calculating the estimated improvement of a group over another (in terms of clicks per user). Much like a A/B-test, but I'm using PyMC to be nice and Bayesian about it.

This data and code works like a charm.

import pymc as pm
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

### Here's the data
# Bins of users that clicked 1 time, 2 times, 3 times, ...
bins_A = [207763, 223077, 210613, 181571, 168385, 159171, 146068, 128502,
110505,  94379,  79315,  67084,  57527,  48867,  41862]
bins_B = [219812, 228003, 208490, 182409, 173357, 164470, 151033, 132412,
113750,  95835,  81206,  67876,  58057,  49005,  41808]

clicks = range(1, len(bins_A))

# Start with uniform probability over the bins
p_A = pm.Dirichlet("p_A", theta=np.ones(len(bins_A)))
p_B = pm.Dirichlet("p_B", theta=np.ones(len(bins_B)))

# A multimodal dist. using the probabilitys of bins
obs_A = pm.Multinomial("obs_A", p=p_A, n=sum(bins_A), value=bins_A, observed=True)
obs_B = pm.Multinomial("obs_B", p=p_B, n=sum(bins_B), value=bins_B, observed=True)

@pm.deterministic
def percent_better(p_B=p_B, p_A=p_A, clicks=clicks):
exp_clicks_B = np.dot(p_B.astype(float)/sum(p_B), clicks)
exp_clicks_A = np.dot(p_A.astype(float)/sum(p_A), clicks)

return ((exp_clicks_B / exp_clicks_A) - 1)*100.0

model = pm.Model([p_A, p_B,
obs_A, obs_B,
percent_better])

map_ = pm.MAP(model)
map_.fit()
mcmc = pm.MCMC(model)
mcmc.sample(35000, burn=25000, thin=2)

percent_better_samples = mcmc.trace("percent_better")[:]

print "Probability B > A: {}".format((percent_better_samples > 0).mean())
print "Confidence interval of B:s lift over A:"
print np.percentile(percent_better_samples, 2.5)
print np.percentile(percent_better_samples, 97.5)

print "MCMC error: {}".format(mcmc.stats()['percent_better']['mc error'])
pm.Matplot.plot(mcmc)


## The good results

Nice smooth distribution of possible values of how much group B is better than group A:

## The wierd results

But running it with this data:

bins_A = [1750102,  286721,  122232,   53109,   35203,   23628,   16135,
18991,   24309,   11363,    9732,    8494,    5911,    4374,
3526,    2462,    2186,    1909,    1811,    1684]
bins_B = [1726921,  279424,  111627,   48393,   29513,   20356,   13086,
18364,   23361,   10805,    8752,   10323,    6007,    4252,
3039,    2172,    1829,    1670,    1617,    1569]


There's no convergence and the calculations seems to blow up:

I tried all sorts of settings and read Bayesian Methods for Hackers carefully without luck. Any ideas would be much appreciated!

• I would assume that the number of users who clicked $n$ times would decrease as $n$ grow larger, thus modeling the distribution as $exp$ makes more sense – Uri Goren Mar 24 '16 at 19:26
• @UriGoren Wouldn't that perhaps make it more sensitive to outliers? In this previous question you can se how the data is distributed: http://stats.stackexchange.com/questions/187051/... exp or even lognorm does not really fit that well. That's why I settled for multinomal. – cowboyvspirate Mar 24 '16 at 20:36
• @UriGoren Now you got me thinking! Currently my prior is same probability for every bin in the multinomal. This should clearly be set to something like a exp to be a much more informative prior. I'll try that! – cowboyvspirate Mar 24 '16 at 20:42
• Also, since Dirrechlet is the conjugate of the Multinomial, why do you need pymc for ? – Uri Goren Mar 24 '16 at 20:52
• @UriGoren You're right, but I like it as a framework and I tried lots of different models for this. – cowboyvspirate Mar 24 '16 at 21:09