I am building a bayesian A/B testing code in python to compare conversion rates between two samples which I will call $A$ and $B$. My aim is to produce a mean conversion rate, in percentage, for $A$ and a mean conversion rate, in percentage, for $B$ then report a gain of $A$ over $B$ with an error on the gain.

The code is pretty straightforward and this is not a coding question but I will poste it anyway:

from scipy.stats import beta

# class for bayesian AB testing
class BayesABTest:
    # get beta distribution function from total players and those who came back
    def getBetaDistribution(self, total, cameBack):
        a = cameBack + 1
        b = total- cameBack + 1
        return beta(a, b)
    # calculate gain between two samples for active and passive players
    def getGain(self, totalA, cameBackA, totalB, cameBackB):
        beta1 = self.getBetaDistribution(totalA, cameBackA)
        beta2 = self.getBetaDistribution(totalB, cameBackB)
        gain = (beta1.mean() - beta2.mean())*100
        return gain, beta1, beta2

Now the gain in percentage is pretty simple to compute as the difference between the means of the two beta distributions. The error, however, is a bit less obvious to me. In a usual setting, I would compute the two standard errors on the means of each sample, add them, then take the square root like so:

$$ SE_{diff} = \sqrt{\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}} $$

But since I transformed my means into distributions, I am not sure how this works. So my question are:

  • How can I compute the standard error on a beta distribution?
  • Can I simply report the error as the square root of the sum of the two standard deviations of the beta distributions squared like so : $error = \sqrt{s_1^2 + s_2^2}$

Thanks for the help!

  • 1
    $\begingroup$ Here's a computational way to get what you want: You have the posterior distribution of the sucess probability p (which is beta(a, b) for each variant A, B. Simulate random values from both posterior distributions (say 10000 each), getting two arrays of 10,000 samples each for A and B. Then, take the difference of these two arrays to get the posterior distribution for pr(success A) - pr(success B), from which you can take the empirical mean, standard deviation, or any other quantity of interest you'd like. $\endgroup$
    – aranglol
    Nov 4, 2022 at 17:21
  • $\begingroup$ You could just simulate it. Since you're using a conjugate prior, just sample from each posterior and estimate the standard error through sampling. $\endgroup$ Nov 4, 2022 at 17:36
  • $\begingroup$ Thank you both for your input. @aranglol I tried simulating random values from both posterior distributions then compute the difference as you suggest. It seems to me that the parameters of this new posterior distribution for pr(success A) - pr(success B) depend entirely on the number of random values simulated... I can get either a very low standard error by doing so or a very large one if I only take, say, 100 random values. Any thoughts on this? $\endgroup$
    – bmasri
    Nov 7, 2022 at 9:29
  • $\begingroup$ Also, can is it not safe to assume that a standard deviation in a Bayesian posterior such as mine is equal to the standard error? Such thing would imply that I can simply consider that simulating n random variables from both posteriors and taking the difference then taking the standard deviation will suffice to report an error on the difference in means. $\endgroup$
    – bmasri
    Nov 7, 2022 at 9:58
  • $\begingroup$ @bmasri What you are observing is Monte Carlo error, not the error in your parameter that you are trying to estimate. You should start to see convergence in the estimated standard deviation for an arbitrarily large number of samples. For your second question, yes! I think there might be some confusion with frequentist vs. Bayesian statistics here. In Bayes land, there isn't really a notion of "standard error" because that speaks to a sampling distribution, a frequentist idea. We don't work with sampling distributions, because you have the exact distribution already based on your model and... $\endgroup$
    – aranglol
    Nov 7, 2022 at 17:34

1 Answer 1


I will outline my comment as an answer.

First, some clarification. In Bayesian statistics, there is no notion of "standard error" in the context of trying to conduct inference on a statistic. The standard error of an estimator is the standard deviation of the estimator's sampling distribution. A sampling distribution relies on the notion of repeated random sampling - a frequentist interpretation of probability.

Rather, what you want is the posterior distribution of the differences in parameters - specifically the posterior distribution of $W := \text{P(success A)} - \text{P(success B)}$.

Once you have the distribution of W, you can describe its qualities in whatever capacity you wish (credible intervals, i.e. quantiles, highest density intervals, etc.) and you can get point estimates by reporting the mean, median, etc.

In theory, you can derive W by thinking about the distribution of W = X - Y, where X and Y are independent Beta random variables. This might give you some distribution in closed form, if you can work out the math.

A much easier way to do this is described in the comments - use Monte Carlo simulation. If you randomly sample an arbitrarily large number of realizations from each posterior/beta distribution (say 10000 each), and simply take the difference between these two arrays, you will have obtained 10000 realizations of W! From here, you can get a Monte Carlo estimate of whatever quantity you wish to conduct inference on W - mean, median, ... using empirical estimators.

In the comments, you stated concerns regarding the "standard error decreasing for a large nunber of random realizations". This is Monte Carlo error - the error that you obtain by approximating some statistic of W with samples from W (as opposed to using math to work out the distribution and/or quantities of W directly). This is error that is unrelated to what is implied by your data or prior beliefs, but rather the "cost" of using simulation. Luckily, as you already noticed, the error goes to 0 for a large number of random samples, so as long as you simulate a large number of realizations, this error will be negligible.

To see this: take the mean of W, derived via. simulation, and compare it to the mean difference that you have already calculated in your code. What you should find is that your Monte Carlo estimate for the mean of W approaches the estimate that you explicitly calculated in closed form for a progressively larger number of samples.

  • 1
    $\begingroup$ +1 for the clear answer and explanations. Indeed, I have observed that with increasing number of random variables drawn from each distribution $X$ and $Y$, the empirically calculated parameters of the distribution of $W=X-Y$ are the same as the ones calculated explicitly in closed form using scipy's methods. Thanks again! $\endgroup$
    – bmasri
    Nov 8, 2022 at 9:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.