Note: Thinking again, if you use the normal approximation (as I did below), deriving the mass function as done below is a bit of overkill - all you need are the first two moments of $Q = \sum_{i} Z_{i}$, which, by properties of means and variances are simply $nE(X_{i}-Y_{i})$ and $n{\rm var}(X_{i}-Y_{i})$. The mass function would be unecessary if you wanted to do an exact test, which would require the CDF of $Q$, though.
So, you're basically asking about the distribution of $Q = \sum_{i} Z_{i}$ where $Z_{i} = Y_{i} - X_{i}$ so you can form a confidence interval for the expected value of $Q$, $E(Q)$, where $Y_{i} \sim {\rm Binomial}(n_{1}, p_{1})$ and $X_{i} \sim {\rm Binomial}(n_{2}, p_{2})$ for each $i$, right? It wasn't clear to me whether or not $n_{1} = n_{2}$, so I will not assume this is true.
From the law of total probability, we have, for each $i$,
$$ P(X_{i} - Y_{i} = k) = E_{Y} \Big( P(X_{i}-Y_{i} = k) \Big) = \sum_{y=0}^{\infty} P(Y_{i}=y)P(X_{i} = k+y) $$
so
$$ P(X_{i} - Y_{i} = k) = \sum_{y=0}^{n_{1}}
{n_{1} \choose y} p_{1}^{y} (1-p_{1})^{n_{1} - y}
{n_{2} \choose k+y} p_{2}^{k+y} (1-p_{2})^{n_{2} - k - y}
$$
Where ${n_{2} \choose k+y}$ is defined as 0 when $k+y > n_{2}$. Using this formula, one can easy, numerically, calculate the mean and variance of $X_{i} - Y_{i}$ using, for example, this R code:
massfcn <- function(k,n1,n2,p1,p2)
{
s <- 0
for(y in 0:n1)
{
term1 <- choose(n1,y)*(p1^y)*((1-p1)^(n1-y))
term2 <- choose(n2,k+y)*(p2^(k+y))*((1-p2)^(n2-y-k))
s <- s + term1*term2
}
return(s)
}
# expected value of X-Y for given n1,n2,p1,p2
E <- function(n1,n2,p1,p2)
{
# range of possible values for X-Y
s = seq(-n1,n2,by=1)
len = length(s)
mn = 0
for(i in 1:len) mn = mn + massfcn(s[i],n1,n2,p1,p2)*s[i]
return(mn)
}
# variance of X-Y for given n1,n2,p1,p2
VAR <- function(n1,n2,p1,p2)
{
# range of possible values for X-Y
s = seq(-n1,n2,by=1)
len = length(s)
mn = 0
for(i in 1:len) mn = mn + massfcn(s[i],n1,n2,p1,p2)*(s[i]^2)
return(mn - E(n1,n2,p1,p2)^2)
}
You can use the functions above to obtain $\mu = E(X_{i} - Y_{i})$ and $\sigma^{2} = {\rm var}(X_{i} - Y_{i})$. Then if each of the $i$'s are independent, it follows that $E(Q) = n \mu$ and ${\rm var}(Q) = n \sigma^{2}$. What you're looking to test is equivalent to testing whether $\mu = 0$ - scaling by $n$ to produces the sample mean, $ \overline{Z} = \frac{1}{n} Q$, which has mean $\mu$ and variance $\sigma^{2}/n$. Therefore, rejecting $\mu = 0$ if the interval
$$ (\overline{Z} - 1.96\sqrt{\sigma^{2}/n}, \overline{Z} + 1.96\sqrt{\sigma^{2}/n}) $$
contains 0 is an approximate $5\%$ test of that null hypothesis (approximate in the sense that the central limit theorem is used to form this interval). Technically, since $p_{1}, p_{2}$ are not known and are estimated from the data, the $\sigma^{2}$ in the above equation should be replaced by $\hat{\sigma}^{2}$, but the procedures and the large sample distribution remains the same.
# example
# some (n=50) fake binomal(12,.3) and binomial(12,.4) data.
x = rbinom(50, 12, .3)
y = rbinom(50, 12, .4)
# numbers of trials
n1 = 12
n2 = 12
# sample proportions
p1 = mean(y)/n1
p2 = mean(x)/n2
# estimate E(X-Y)
MEAN = E(n1,n2,p1,p2)
# estiamte var(X-Y)
VARIANCE = VAR(n1,n2,p1,p2)/50
# get an approximate confidence interval
CI = c(MEAN-1.96*sqrt(VARIANCE), MEAN+1.96*sqrt(VARIANCE))
