How can we incorporating uncertainty about our data into Bayesian inference? I want to use a Bayesian approach to estimate the parameter $\theta$ of a binomial distribution $\textsf{Binomial}(\theta,n)$ with the number of Bernoulli experiments $x_i \in \{0,1\}$ being known and fixed to $n$.
In the textbook example, the binomial likelihood only takes the sum of positive vs. negative outcomes of $x_i$ into account:
\begin{equation}
f(x|\theta) = \theta^{\sum x_i} (1-\theta)^{n-\sum x_i}
\end{equation}
Together with a uniform prior $f(\theta)=\textsf{Uniform}(0,1)$, the posterior becomes the Beta distribution
\begin{equation}
f(\theta|s) = \textsf{Beta}(\theta|s+1,n-s+1) . 
\end{equation}
, where $s=\sum x_i$. Unfortunately, $x_i$ is not directly observable. Instead, we are given $\pi_i$ as features to e.g. a neural network, which then provides a probabilistic mapping $f(x_i|\pi_i)$. This probability represents uncertainty about the true, latent value $x_i$.
How can we incorporate this uncertainty about the data into the estimate? E.g. if $f(x_i|\pi_i)=0.51\ \forall\ 0\leq\ i<n$, we would expect the MAP estimate of $\theta$ to be a lot further off from $1$ than if $f(x_i|\pi_i)=0.91$ throughout.
 A: The distribution for $\theta$ conditional on $s$ (and $n$) is
\begin{equation}
p(\theta|s) = \textsf{Beta}(\theta|s+1,n-s+1) . 
\end{equation}
The problem is that we don't observe $s$.
Instead we are given $\pi = (\pi_1, \ldots, \pi_n)$, where
\begin{equation}
p(x_i|\pi_i) = \textsf{Bernoulli}(x_i|\pi_i)  
\end{equation}
and
\begin{equation}
p(x|\pi) = \prod_{i=1}^n p(x_i|\pi_i) .
\end{equation}
We can compute $p(s|\pi)$ from $p(x|\pi)$ and then
\begin{equation}
p(\theta|\pi) = \sum_{s=0}^n p(\theta|s)\,p(s|\pi) ,
\end{equation}
which is a mixture of beta distributions.
The problem now is that $p(s|\pi)$ is the Poisson-binomial distribution which is difficult to compute analytically unless $n$ is quite small. One solution is to use simulation to approximate the distribution. In particular, let
\begin{equation}
p(\theta|\pi) \approx \frac{1}{R} \sum_{r=1}^R p(\theta|s^{(r)}) ,
\end{equation}
where $s^{(r)} \sim p(s|\pi)$.
Let me provide more detail. Let
\begin{equation}
s^{(r)} = \sum_{i=1}^n x_i^{(r)} ,
\end{equation}
where
\begin{equation}
x_i^{(r)} \sim p(x_i|\pi_i) .
\end{equation}
Now define the approximate mixture weights:
\begin{equation}
\hat w_j = \frac{1}{R}\sum_{r=1}^R 1(s^{(r)} = j) ,
\end{equation}
where the indicator function is
\begin{equation}
1(x) = \begin{cases}
1 & \text{$x$ is true} \\
0 & \text{$x$ is false}
\end{cases} .
\end{equation}
Then
\begin{equation}
p(\theta|\pi) \approx \sum_{j=0}^n \hat w_j\,p(\theta|s=j) . 
\end{equation}
You will need to experiment to see how large $R$ needs to be so that you get a sufficiently accurate answer. Depending on $n$ and $\pi$, you might need $R = 10^5$ or more.
