Approximation of product of Bernoulli with different proportions

Question

I want to update a variable $Y$ with Beta (uniform for simplicity, $Y \sim U(0, 1)$) distribution, with Bernoulli information each period... But each period the proportion parameter of the Bernoulli varies. I have a product of Bernoulli densities with different proportions $p_i$. (success $x_i = 1$ or failure $x_i = 0$).

It looks like this:

$ f(y|x_i) \propto \Pi_{i=1}^{n} (p_i y)^{x_i}(1-p_i y)^{1-x_i} = (p_1 y)^{x_1} (1-p_1 y)^{1-x_1} ... \cdot \ (p_n y)^{x_n} (1-p_n y)^{1-x_n}$
where $y \in [0, 1]$

I would like to approximate it by: $(Py)^{\sum_{i=1}^{n} x_i} (1-Py)^{n - \sum_{i=1}^{n} x_i}$ where $P$ would be a function of $(p_1, p_2, ..., p_n, \sum_{i=1}^{n} x_i, n)$.
I need such an approximation because I would need to summarize the information in all the $p_i$ in a small number of parameters. If needed, I can use at most two parameters, something like:
$(P_1 y)^{\sum_{i=1}^{n} x_i} (1-P_2y)^{n - \sum_{i=1}^{n} x_i}$ where $P_1$ and $P_2$ would be a function of $(p_1, p_2, ..., p_n, \sum_{i=1}^{n} x_i, n)$.

I do not have a good math background, is there any good analytical approximation for this that I am not aware of? Or do I have to solve it numerically ?

---

Edit: Partial solution.

As suggested by Jim in the comment below, I simply have the two closed forms:

$P_1 = (\Pi^n_{i=1} p_i^{x_i})^{1/\sum x_i}$
$P_2 = \frac{1}{y}(1 - (\Pi^n_{i=1} (1 - y p_i)^{1-x_i})^{1/(n - \sum x_i)})$

From this, I can build recursive formulation: from a given $P_1$, knowing the present $x_i$ and $p_i$, I can build the next $P_1$.

However, it does not work as well for $P_2$. As you can see, $P_2$ depends on $y$. So I actually have a distribution of $P_2$... While I would like to gather all the information about $P_2$ in ONE value. This is not feasible (at least I did not find any solution). However, $P_2 \times y$ seems to be a quasi-linear function of $y$: I think I will save the coefficient of the linear regression as my coefficient of information $P_2$ (by default). This is the best approximation I can think of, except if any of you have a better idea.

Answer 1

If your endgame is to get the posterior density for y, I'm not understanding why you need an approximation when you can get the exact posterior density for any particular set of data. Here's how to do it using Mathematica:

(* Data *)
n = 10;
p = {2/10, 3/10, 222/1000, 15/100, 19/100, 4/10, 44/100, 21/100, 34/100, 1/10};
x = {1, 1, 0, 0, 1, 1, 1, 1, 0, 1};

(* Constant of integration *)
const = Integrate[Expand[Product[y^x[[i]] p[[i]]^x[[i]] (1 - y p[[i]])^(1 - x[[i]]), {i, n}]], {y, 0, 1}]
(* 200262223/78125000000000 *)

(* Posterior density for y *)
f[y_] := (1/const) Product[y^x[[i]] p[[i]]^x[[i]] (1 - y p[[i]])^(1 - x[[i]]), {i, n}]

(* Plot density *)
Plot[f[y], {y, 0, 1}, PlotRange -> All, Frame -> True, 
 FrameLabel -> {{Style["Density", Large], ""}, {Style["y", Large],""}}]

(* Check on density integrating to 1 *)
Integrate[f[y], {y, 0, 1}]
(* 1 *)