Bayesian A/B Testing of Count Data

By Evan Miller

June 8, 2014

As a follow-up to A Formula for Bayesian A/B Testing, here's a formula to predict the result of a test with count data (for example, if you're comparing the number of sales per salesman, or number of sales week over week):

\[ {\rm Pr}(\lambda_1 > \lambda_2) = \sum_{k=0}^{\alpha_1-1} \frac{ (\beta_1 + \beta_2)^{-(k+\alpha_2)} \beta_1^k \beta_2^{\alpha_2} }{(k+\alpha_2)B(k+1,\alpha_2)} \]

Here the \(\alpha\) values represent the total event counts in each group, and the \(\beta\) values represent the exposure, that is, the relative opportunity for events to occur. For example, if the first group was “exposed” to events for twice as long as the second group, you would set \(\beta_1 = 2\) and \(\beta_2 = 1\). (Alternatively, you could set \(\beta_1 = 20\) and \(\beta_2 = 10\); the math works out the same.) In the above equation, \(B\) is the beta function.

A derivation and implementation follow; they both closely mirror the binary-outcome case.


Not for the timid! You may want to skip right to the implementation section.

Let \(\lambda_1\) and \(\lambda_2\) be the Poisson parameter for each group. With a gamma-distributed prior belief, the posterior beliefs are given by:

\[ \begin{array}{lr} \lambda_1 \sim {\rm Gamma}(\alpha_1, \beta_1) \\ \lambda_2 \sim {\rm Gamma}(\alpha_2, \beta_2) \end{array} \]

Using the pdf of the gamma distribution, we can get the total probability that \(\lambda_1\) is greater than \(\lambda_2\) by integrating the joint distribution over all values for which \(\lambda_1 > \lambda_2\):

\[ \begin{equation} {\rm Pr}(\lambda_1 > \lambda_2) = \int_0^\infty \int_{\lambda_2}^\infty \frac{\beta_1^{\alpha_1} \lambda_1^{\alpha_1 - 1} e^{-\beta_1 \lambda_1}}{\Gamma(\alpha_1)} \frac{\beta_2^{\alpha_2} \lambda_2^{\alpha_2 - 1} e^{-\beta_2 \lambda_2}}{\Gamma(\alpha_2)} d\lambda_1 d\lambda_2 \end{equation} \]

Evaluating the inner integral, the equation becomes:

\[ \begin{equation} \label{pr_eval_inner} \int_0^\infty Q(\alpha_1, \beta_1 \lambda_2) \frac{\beta_2^{\alpha_2} \lambda_2^{\alpha_2-1} e^{-\beta_2 \lambda_2}}{\Gamma(\alpha_2)} d\lambda_2 \end{equation} \]

Where \(Q\) is the upper incomplete regularized gamma function. Using the identity \(Q(1,z) = e^{-z}\) and the recursive relationship

\[ \begin{equation} Q(a + n, z) = Q(a, z) + z^a e^{-z} \sum_{k=0}^{n-1} \frac{z^k}{\Gamma(a+k+1)} \end{equation} \]

we can express \(Q\) as:

\[ \begin{equation} \label{igamma_sum} Q(n, z) = e^{-z}\sum_{k=0}^{n-1}\frac{z^k}{\Gamma(k+1)} \end{equation} \]

The probability integral \(\eqref{pr_eval_inner}\) can therefore be written:

\[ \begin{array}{ll} {\rm Pr}(\lambda_1 > \lambda_2) &=& \int_0^\infty e^{-\beta_1 \lambda_2} \left( \sum_{k=0}^{\alpha_1-1} \frac{(\beta_1 \lambda_2)^k}{\Gamma(k+1)}\right) \frac{\beta_2^{\alpha_2} \lambda_2^{\alpha_2-1} e^{-\beta_2 \lambda_2}}{\Gamma(\alpha_2)} d\lambda_2 \\ &=& \sum_{k=0}^{\alpha_1-1} \int_0^\infty \frac{e^{-\beta_1 \lambda_2} (\beta_1 \lambda_2)^k}{\Gamma(k+1)} \frac{\beta_2^{\alpha_2} \lambda_2^{\alpha_2-1} e^{-\beta_2 \lambda_2}}{\Gamma(\alpha_2)} d\lambda_2 \\ &=& \sum_{k=0}^{\alpha_1-1} \frac{\beta_1^k \beta_2^{\alpha_2}}{\Gamma(k+1) \Gamma(\alpha_2)} \int_0^\infty e^{-(\beta_1+\beta_2) \lambda_2} \lambda_2^{k+\alpha_2-1} d\lambda_2 \\ &=& \sum_{k=0}^{\alpha_1-1} \frac{\beta_1^k \beta_2^{\alpha_2}}{\Gamma(k+1) \Gamma(\alpha_2)} (\beta_1 + \beta_2)^{-(k+\alpha_2)} \Gamma(k+\alpha_2) \\ &=& \sum_{k=0}^{\alpha_1-1} \beta_1^k \beta_2^{\alpha_2} (\beta_1 + \beta_2)^{-(k+\alpha_2)} \frac{\Gamma(k+\alpha_2)}{\Gamma(k+1) \Gamma(\alpha_2)} \frac{k+\alpha_2}{k+\alpha_2} \end{array} \]

Using \(\Gamma(z+1) = z\Gamma(z)\), we have:

\[ \begin{equation} \label{gamma_version} {\rm Pr}(\lambda_1 > \lambda_2) = \sum_{k=0}^{\alpha_1-1} \beta_1^k \beta_2^{\alpha_2} (\beta_1 + \beta_2)^{-(k+\alpha_2)} \frac{\Gamma(k+\alpha_2+1)}{\Gamma(k+1) \Gamma(\alpha_2)} \frac{1}{k+\alpha_2} \end{equation} \]

Replacing the gamma functions with a beta function, we have:

\[ \begin{equation} \label{pr_alpha_b} {\rm Pr}(\lambda_1 > \lambda_2) = \sum_{k=0}^{\alpha_1-1} \frac{ (\beta_1 + \beta_2)^{-(k+\alpha_2)} \beta_1^k \beta_2^{\alpha_2} }{(k+\alpha_2)B(k+1,\alpha_2)} \end{equation} \]


Here's a quick implementation of \(\eqref{pr_alpha_b}\) in Julia:

function probability_1_beats_2(α_1, β_1, α_2, β_2)
    total = 0.0
    for k = 0:(α_1-1)
        total += exp(k * log(β_1) + α_2 * log(β_2) - (k+α_2) * log(β_1 + β_2)
            - log(k + α_2) - lbeta(k + 1, α_2))
    return total

The beta function produces very large numbers, so if you’re getting infinite values in your program, be sure to work with logarithms, as in the code above. Your standard library’s log-beta function will come in handy here.

If you don’t have log-beta available, it's easy enough to define one with the log-gamma function and the identity:

\[ \log(B(a, b)) = \log(\Gamma(a)) + \log(\Gamma(b)) - \log(\Gamma(a+b)) \]

You're reading, a random collection of math, tech, and musings. If you enjoyed this, you might also like:

If you run A/B tests frequently, you should check out Evan's Awesome A/B Tools:

Sample Size Calculator

Chi-Squared Test

Two-Sample T-Test

Finally, if you own a Mac, my desktop statistics software Wizard can help you analyze more data in less time and communicate discoveries visually without spending days struggling with pointless command syntax. Check it out!

Statistical analyzer

Back to Evan Miller's home pageFollow on TwitterSubscribe to RSS