# Mismatch between analytical and simulation way for Bayesian estimation of binomial data

In short, Below, I'm asking why the analytical way of getting the posterior for the difference between two proportions (binomial data) and simulation way of the same don't match?

(Note: All code is R code)

# Details

Suppose I have the following information for two groups (g1 and g2) of binomially distributed data:

g1.trials = g2.trials = 100

g1.success = 55                   ;    g2.success = 45
g1.fail = g1.trials - g1.success  ;    g2.fail = g2.trials - g2.success


with prior information for a beta distribution:

alpha1 = alpha2 = 1.1
beta1 =  beta2 = 1.1


My goal is to obtain the posterior probability distribution for $\Delta~(p_1 -p_2)$ i.e., the difference between "true proportion in group 1 - true proportion in group 2".

Simulation solution:

It is widely known the the posterior for each group is:

set.seed(0)
post.p1 = rbeta(1e6, alpha1 + g1.success, beta1 + g1.fail)
post.p2 = rbeta(1e6, alpha2 + g2.success, beta2 + g2.fail)


And their difference i.e., $\Delta~(p_1 -p_2)$ is:

delta.p = post.p1 - post.p2  # the difference

hist(delta.p, freq = FALSE, xlab = bquote(Delta[~(p[1]-p[2])]))

med = median(delta.p)                    # 0.09822197
CI = quantile(delta.p, c(0.025, 0.975))  # -0.03834067  0.23272198


Analytical solution:

We follow the usual Bayesian way (i.e., for each group: prior x likelihood):

prior1 = function(x) dbeta(x, alpha1, beta1);
like1 = function(x) dbinom(g1.success, g1.trials, x);
k1 = integrate(function(x) prior1(x)*like1(x), 0, 1)[[1]];
post1 = function(x) prior1(x)*like1(x) / k1;

prior2 = function(x) dbeta(x, alpha2, beta2);
like2 = function(x) dbinom(g2.success, g2.trials, x);
k2 = integrate(function(x) prior2(x)*like2(x), 0, 1)[[1]];
post2 = function(x) prior2(x)*like2(x) / k2;

kdelta = integrate(function(x) post1(x) - post2(x), 0, 1)[[1]];

postdelta = function(x) post1(x) - post2(x) / kdelta;

curve(postdelta, xlab = bquote(bold(Delta[~(p[1]-p[2])]))) # completely different!
(mode = optimize(postdelta, 0:1, maximum = TRUE)[[1]]) # 0.4500863 Completely different
# from simulation median!!!!!


The distribution of the difference of two random variables $X$ and $Y$, $X-Y$, is not the difference of their distributions, $f_X-f_Y$.
The general derivation of the density of $X-Y$ is the convolution formula $$f_{X-y}(z)=\int_\mathcal{X}f_X(x)f_y(x-z)\text{d}x$$
In the case of two Beta variates, $X\sim\mathcal{B}e(\alpha+n_1,\beta+n_2)$ and $Y\sim\mathcal{B}e(\alpha+m_1,\beta+m_2)$ the density of the difference $X-Y$ is thus provided by \begin{align*}f_{X-y}(z)&=\int_\mathcal{X}f_X(x)f_y(x-z)\text{d}x\\ &=\text{B}(\alpha+n_1,\beta+n_2) \text{B}(\alpha+m_1,\beta+m_2)\\ &\times\int_{\max\{0,z\}}^{\min\{1,1+z\}} x^{\alpha+n_1-1}(1-x)^{\beta+n_2-1} (x-z)^{\alpha+m_1-1}(1-x+z)^{\beta+m_2-1}\text{d}x\end{align*} for which there is a closed form expression when $\alpha$ and $\beta$ are integers, since the integrand is then a polynomial in $x$.