Reporting standard error on gain using Bayesian A/B testing 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 and 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!
 A: 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.
