I'm working on a question with simultaneous learning about an underlying population and individual members of the population. The basic setup is:
Let $N_g$ be the size of a population.
At any point in time $t$, each object $i \in \{1, ... , N_g\}$ has been sampled $\{n_1, ..., n_{N_g}\}$ times.
For each object $i$, we observe $r_i$ successful realizations out of the $n_i$ observations. With this in mind, I'm using the Beta-binomial distribution so that: \begin{align*} r_i &\sim Bin(n_i, \theta_i) \\ \theta_i &\sim Beta(\alpha,\beta) \\ (\alpha,\beta) &\sim p_0(\alpha,\beta) \end{align*}
Given a sequence of observations $\{r_i,n_i\}_{i=1}^{N_g}$, I'm hoping to calculate the marginal posterior expectations: \begin{align*} \mathbb{E}(\theta_i \vert \{r_i,n_i\}_{i=1}^{N_g}) \end{align*}
That is, I'm hoping to calculate the expected probability of success for object $i$ given the entire sample from objects $1$ to $N_g$.
My logic in proceeding has been to calculate: \begin{align*} Pr(\{\theta_{j},r_j,n_j\}_{j^{\prime}}, \alpha,\beta) &= \left[\prod_{j=1}^{N_{g}} {n_{j}\choose r_{j}} \theta_{j}^{r_{j}}(1-\theta_{j})^{n_{j}-r_{j}}\right] \cdot \left[\prod_{j=1}^{N_{g}} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta_{j}^{\alpha-1}(1-\theta_{j})^{\beta-1} \right] \cdot p^{g}_{i}(\alpha,\beta) \\ \text{given a uniform prior } p_0(\alpha,\beta): \\ Pr(\alpha,\beta \vert \{r_j,n_j\}_{j^{\prime}}) & \propto \left[\prod_{j=1}^{N_{g}} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \frac{\Gamma(\alpha + r_{j})\Gamma(\beta+n_{j}-r_{j})}{\Gamma(\alpha+\beta + n_{j})} \right] \\ \text{so:} \\ Pr(\theta_{j} \vert \{r_j,n_j\}_{j^{\prime}}) &= \int_{(\alpha,\beta)} Pr(\theta_{j} \vert \alpha,\beta\{r_j,n_j\}_{j^{\prime}}) \cdot Pr(\alpha,\beta \vert \{r_j,n_j\}_{j^{\prime}}) d(\alpha,\beta) \\ \Rightarrow \mathbb{E}(\theta_{j} \vert \{r_j,n_j\}_{j^{\prime}}) &= \int_{\alpha,\beta} \left[\frac{\alpha + r_j}{\alpha+\beta+n_j}\right] \cdot \prod_{j=1}^{N_{g}} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \frac{\Gamma(\alpha + r_{j})\Gamma(\beta+n_{j}-r_{j})}{\Gamma(\alpha+\beta+n_{j})} \cdot Cd(\alpha,\beta) \end{align*}
Unfortunately, I'm getting stuck in this calculation. I've been trying to find literature where someone has already worked this out (I do Economics, and we're often 40+ years behind math, so I'm sure it's been done.), but I'm coming up empty.
I've been trying to read through discussions of the Dirichlet distribution (insofar as the beta-binomial is a subcase of Dirichlet), but I can't seem to be able to find anything about this type of posterior.
Does anyone have any idea / recommendation where I could look?
Thank you!