I've got a series of probabilities estimated as the fitted values of a logistic regression. Call them $\hat y_i$. I've got them on the logit scale, so one can get the estimated probability as $(1+e^{-\hat y_i})^{-1}$. The logit model is somewhat uncertain; I've also got standard errors of $\hat y_i$, call them $se(\hat y_i)$.
Given both uncertainty in the outcome of the Bernouli trials, and uncertainty in each Bernouli trial's probability, what is the probability mass function for the sum of the many Bernouli trials? Is this possible to get it analytically?
To clarify: I start with a logistic regression model $$ pr(y = 1) = logit(\mathbf{X'\beta}+\epsilon) $$ which gets me $\mathbf{X'\hat\beta} = \hat y$ on the logit scale. Each $\hat y_i$ is associated with different $\mathbf{X}_i$, so will yield a different estimated probability $\widehat{pr}_i = (1+e^{-\hat y_i})^{-1}$. In addition, each $\hat y$ is associated with a standard error, so that each $\hat y_i$ is (asymptoticaly) $\mathcal{N}(\hat y, se(\hat y))$ on the logit scale. So the confidence intervals on the probabilities aren't symmetrical unless $\hat y = 0$.
OK. So you see that there is uncertainty about whether $pr_i$ will manifest as a one or a zero. You also see that there is uncertainty about the $pr_i$ itself. Now, say I have $i = 1,2,...,n$. What is the distribution of $$ \displaystyle\sum_n y_i $$ which is an integer between zero and $n$? It will be the sum of many bernouli trials, each of which has a different probability, of which each has a different distribution.
Edit: Sorry for any confusion re: the writing of $\hat y$ as being on the logit scale and $y$ as being the $0/1$ outcome.