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)


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:

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$.

  • $\begingroup$ Xi'an, OP seems to wanted to know more that what you have very briefly provided. It will be useful to provide a fuller response focused on Beta random variables. $\endgroup$ – Reza Dec 12 '17 at 22:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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