Consider two independent discrete random variables $y_1$ and $y_2$, both distributed with a Beta-Binomial distribution, with different number of successes $n_1$ and $n_2$ but the same parameters $a$ and $b$

$ p(y_1|n_1,a,b) = {n_1 \choose y_1} \dfrac{B(y_1 + a,n_1 -y_1 +b)}{B(a,b)} $

$ p(y_2|n_2,a,b) = {n_2 \choose y_2} \dfrac{B(y_2 + a,n_2 -y_2 +b)}{B(a,b)} $

Consider a discrete variable $Z = y_1 + y_2$. Is $Z$ distributed as well as a Beta-Binomial (with parameters $n_1+n_2$, $a'$ and $b'$?

I could not prove it in an analytic form so far, but I have been trying it out with some simulations. Reparametrising $a$ and $b$ as $\mu = \dfrac{a}{a+b}$ and $\rho = \dfrac{1}{a+b+1}$, $y_1$ and $y_2$ have mean $\mu n_1$ and $\mu n_2$ respectively and variance $\mu(1-\mu)n_1(1 + (n_1-1)\rho)$ and $\mu(1-\mu)n_2(1 + (n_2-1)\rho)$ respectively.

Based on independence, the mean of $Z$ is $\mu(n_1+n_2)$ and the variance of $Z$ is $\mu(1-\mu)(n_1(1 + (n_1-1)\rho) + n_2(1 + (n_2-1)\rho))$. If $Z$ was distributed according to a beta-binomial distribution, then it would have paramters $\mu'$ and $\rho'$, with $\mu'= \mu$ and

$\rho' = \dfrac{\dfrac{n_1(1 + (n_1-1)\rho) + n_2(1 + (n_2-1)\rho)}{n_1+n_2} - 1}{n_1+n_2-1} = \rho \dfrac{n_1(n_1-1)+n_2(n_2-1)}{(n_1+n_2)(n_1+n_2 -1)} $

Here is some code to generate $Z$ as a sum of two independent Beta-Binomials (sorry about the code, R is not my main language)

n1 = 20
n2 = 50
mu = .6
k  = 20

p1  = rbeta(1e6,mu*k,(1-mu)*k)
y1  = rbinom(1e6,n1,p1)

p2  = rbeta(1e6,mu*k,(1-mu)*k)
y2  = rbinom(1e6,n2,p2)

z   = y1+y2

rho  = 1/(k+1)
rho1 = rho*(n1*(n1-1)+n2*(n2-1))/((n1+n2)*(n1+n2-1))
k1   = 1/ rho1 - 1
p3  = rbeta(1e6,mu*k1,(1-mu)*k1)
z1  = rbinom(1e6,n1+n2,p3)

print(c(var(z),var(z1)))
plot(density(z,width= 3))
lines(density(z1,width = 3))

Consider two independent discrete random variables $y_1$ and $y_2$, both distributed with a Beta-Binomial distribution, with different number of successes $n_1$ and $n_2$ but the same parameters $a$ and $b$

$ p(y_1|n_1,a,b) = {n_1 \choose y_1} \dfrac{B(y_1 + a,n_1 -y_1 +b)}{B(a,b)} $

$ p(y_2|n_2,a,b) = {n_2 \choose y_2} \dfrac{B(y_2 + a,n_2 -y_2 +b)}{B(a,b)} $

Consider a discrete variable $Z = y_1 + y_2$. Is $Z$ distributed as well as a Beta-Binomial (with parameters $n_1+n_2$, $a'$ and $b'$?

I could not prove it in an analytic form so far, but I have been trying it out with some simulations, at least to check whether the assumption is wrong in some cases. Reparametrising $a$ and $b$ as $\mu = \dfrac{a}{a+b}$ and $\rho = \dfrac{1}{a+b+1}$, $y_1$ and $y_2$ have mean $\mu n_1$ and $\mu n_2$ respectively and variance $\mu(1-\mu)n_1(1 + (n_1-1)\rho)$ and $\mu(1-\mu)n_2(1 + (n_2-1)\rho)$ respectively.

Based on independence, the mean of $Z$ is $\mu(n_1+n_2)$ and the variance of $Z$ is $\mu(1-\mu)(n_1(1 + (n_1-1)\rho) + n_2(1 + (n_2-1)\rho))$. If $Z$ was distributed according to a beta-binomial distribution, then it would have paramters $\mu'$ and $\rho'$, with $\mu'= \mu$ and

$\rho' = \dfrac{\dfrac{n_1(1 + (n_1-1)\rho) + n_2(1 + (n_2-1)\rho)}{n_1+n_2} - 1}{n_1+n_2-1} = \rho \dfrac{n_1(n_1-1)+n_2(n_2-1)}{(n_1+n_2)(n_1+n_2 -1)} $

Here is some code to generate $Z$ as a sum of two independent Beta-Binomials (sorry about the code, R is not my main language)

n1 = 20
n2 = 50
mu = .6
k  = 20

p1  = rbeta(1e6,mu*k,(1-mu)*k)
y1  = rbinom(1e6,n1,p1)

p2  = rbeta(1e6,mu*k,(1-mu)*k)
y2  = rbinom(1e6,n2,p2)

z   = y1+y2

rho  = 1/(k+1)
rho1 = rho*(n1*(n1-1)+n2*(n2-1))/((n1+n2)*(n1+n2-1))
k1   = 1/ rho1 - 1
p3  = rbeta(1e6,mu*k1,(1-mu)*k1)
z1  = rbinom(1e6,n1+n2,p3)

print(c(var(z),var(z1)))
plot(density(z,width= 3))
lines(density(z1,width = 3))

I have been trying this code for different values of $n_1$, $n_2$, $\mu$ and $k$, but in all the cases the variances using the sum of two beta-binomials or an appropriately tuned beta-binomial are very similar (the densities look indistinguishable)

Consider two independent discrete random variables $y_1$ and $y_2$, both distributed with a Beta-Binomial distribution, with different number of successes $n_1$ and $n_2$ but the same parameters $a$ and $b$

$ p(y_1|n_1,a,b) = {n_1 \choose y_1} \dfrac{B(y_1 + a,n_1 -y_1 +b)}{B(a,b)} $

$ p(y_2|n_2,a,b) = {n_2 \choose y_2} \dfrac{B(y_2 + a,n_2 -y_2 +b)}{B(a,b)} $

Consider a discrete variable $Z = y_1 + y_2$. Is $Z$ distributed as well as a Beta-Binomial (with parameters $n_1+n_2$, $a'$ and $b'$?

I could not prove it in an analytic form so far, but I have been trying it out with some simulations. Reparametrising $a$ and $b$ as $u = \dfrac{a}{a+b}$ and $\rho = \dfrac{1}{a+b+1}$, $y_1$ and $y_2$ have mean $\mu n_1$ and $\mu n_2$ respectively and variance $mu(1-mu)n_1(1 + (n1-1)\rho)$ and $mu(1-mu)n_2(1 + (n2-1)\rho)$ respectively.

Based on independence, the mean of $Z$ is $\mu(n_1+n_2)$ and the variance of $Z$ is $mu(1-mu)(n_1(1 + (n1-1)\rho) + n_2(1 + (n2-1)\rho))$. If $Z$ was distributed according to a beta-binomial distribution, then it would have paramters $\mu'$ and $\rho'$, with $\mu'= \mu$ and

$\rho' = \dfrac{\dfrac{n_1(1 + (n1-1)\rho) + n_2(1 + (n2-1)\rho)}{n1+n2} - 1}{n1+n2-1} = \rho \dfrac{n_1(n_1-1)+n_2(n_2-1)}{(n_1+n_2)(n_1+n_2 -1)} $

Here is some code to generate $Z$ as a sum of two independent Beta-Binomials (sorry about the code, R is not my main language)

n1 = 20
n2 = 50
mu = .6
k  = 20

p1  = rbeta(1e6,mu*k,(1-mu)*k)
y1  = rbinom(1e6,n1,p1)

p2  = rbeta(1e6,mu*k,(1-mu)*k)
y2  = rbinom(1e6,n2,p2)

z   = y1+y2

rho  = 1/(k+1)
rho1 = rho*(n1*(n1-1)+n2*(n2-1))/((n1+n2)*(n1+n2-1))
k1   = 1/ rho1 - 1
p3  = rbeta(1e6,mu*k1,(1-mu)*k1)
z1  = rbinom(1e6,n1+n2,p3)

print(c(var(z),var(z1)))
plot(density(z,width= 3))
lines(density(z1,width = 3))

Consider two independent discrete random variables $y_1$ and $y_2$, both distributed with a Beta-Binomial distribution, with different number of successes $n_1$ and $n_2$ but the same parameters $a$ and $b$

$ p(y_1|n_1,a,b) = {n_1 \choose y_1} \dfrac{B(y_1 + a,n_1 -y_1 +b)}{B(a,b)} $

$ p(y_2|n_2,a,b) = {n_2 \choose y_2} \dfrac{B(y_2 + a,n_2 -y_2 +b)}{B(a,b)} $

Consider a discrete variable $Z = y_1 + y_2$. Is $Z$ distributed as well as a Beta-Binomial (with parameters $n_1+n_2$, $a'$ and $b'$?

I could not prove it in an analytic form so far, but I have been trying it out with some simulations. Reparametrising $a$ and $b$ as $\mu = \dfrac{a}{a+b}$ and $\rho = \dfrac{1}{a+b+1}$, $y_1$ and $y_2$ have mean $\mu n_1$ and $\mu n_2$ respectively and variance $\mu(1-\mu)n_1(1 + (n_1-1)\rho)$ and $\mu(1-\mu)n_2(1 + (n_2-1)\rho)$ respectively.

Based on independence, the mean of $Z$ is $\mu(n_1+n_2)$ and the variance of $Z$ is $\mu(1-\mu)(n_1(1 + (n_1-1)\rho) + n_2(1 + (n_2-1)\rho))$. If $Z$ was distributed according to a beta-binomial distribution, then it would have paramters $\mu'$ and $\rho'$, with $\mu'= \mu$ and

$\rho' = \dfrac{\dfrac{n_1(1 + (n_1-1)\rho) + n_2(1 + (n_2-1)\rho)}{n_1+n_2} - 1}{n_1+n_2-1} = \rho \dfrac{n_1(n_1-1)+n_2(n_2-1)}{(n_1+n_2)(n_1+n_2 -1)} $

Here is some code to generate $Z$ as a sum of two independent Beta-Binomials (sorry about the code, R is not my main language)

n1 = 20
n2 = 50
mu = .6
k  = 20

p1  = rbeta(1e6,mu*k,(1-mu)*k)
y1  = rbinom(1e6,n1,p1)

p2  = rbeta(1e6,mu*k,(1-mu)*k)
y2  = rbinom(1e6,n2,p2)

z   = y1+y2

rho  = 1/(k+1)
rho1 = rho*(n1*(n1-1)+n2*(n2-1))/((n1+n2)*(n1+n2-1))
k1   = 1/ rho1 - 1
p3  = rbeta(1e6,mu*k1,(1-mu)*k1)
z1  = rbinom(1e6,n1+n2,p3)

print(c(var(z),var(z1)))
plot(density(z,width= 3))
lines(density(z1,width = 3))