Calculating confidence interval from binomial distribution I have two variables:
a: 12,13,15,10,9,8
b: 15,15,15,15,15,15

variable 'a' is aggregate numbers of choosing a particular answer (0 or 1) for each case (6 cases in this example). 
So for example for the 1st case 12 out of 15 participants chose 1, 13 out of 15 in the second case, 15 out of 15 in the third etc.
variable 'b' predicts behavioral responses for each case. Here it predicts that for each case the majority of participants should select option 1.
To compare 'a' to 'b' I calculated the difference between each aggregate response: 
diff: 3,2,0,5,6,7 which is a total of: 23

I would like to compute a confidence interval for this score.
 A: 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))

