I'm going to use $\theta$ as the unknown probability of success. With this change, \begin{equation} p(x_i|\theta) = \textsf{Bernoulli}(x_i|\theta) = \theta^{x_i}\,(1-\theta)^{1-x_i} \end{equation} and \begin{equation} p(x|\theta) = \theta^s\,(1-\theta)^{n-s} , \end{equation} where $s = \sum_{i=1}^n x_i$ is the number of successes. Given the uniform prior for $\theta$, $p(\theta) = \textsf{Uniform}(\theta|0,1)$, 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 is that $p(s|\pi)$ is the Poisson-binomial distribution which is difficult to compute 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)$.

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.