Updating a mixture of Beta distributions I have the following problem and I would appreciate any help you could offer. The quantity I'm interested in ($p$) is a proportion and let's say that in my whole population it follows a mixture of beta distributions
$$ p \sim \phi \mathrm{Beta}(\alpha_1 , \beta_1) + (1 - \phi) \mathrm{Beta}(\alpha_2, \beta_2)$$ 
Based on some knowledge on the population, I have starting values for $\phi, \alpha_1, \beta_1, \alpha_2$ and $\beta_2$ but I want to use additional information that arrives sequentially to improve on those starting values. The new information is of the form $p_1, p_2, p_3, ..., p_n$ but I have no knowledge from which of the two $\mathrm{Beta}$ mixing distributions these proportions come from. Furthermore, I don't have the usual Bayesian setting for updating a $\mathrm{Beta}$ prior through the number of trials and success; I only have the sequence of $p_i$.
I don't really know where to start. My gut feeling tells me that based on the value of each new $p_i$, each of the mixing distributions will be updated accordingly, depending on whether it is more likely for $p_i$ to have originated from one or the other distribution. But perhaps I'm completely wrong here so any thoughts on this will be more than welcome. 
So how should I use the new $p_i$ to update the parameters $\phi, \alpha_1, \beta_1, \alpha_2$ and $\beta_2$? And to make things even more complicated, what would happen if the mixing distributions were 3 instead of 2?
 A: I'm going to change the notation a bit to make it easier for me to answer the more general question asked at the end of the original post. I will use $y_i$ instead of $p_i$ to denote the observations and I will use $p$ when expressing probability densities.
The setup for the more general question is this:
\begin{equation}
p(y_i|\phi,\theta) = \sum_{j=1}^k \phi_j\,p(y_i|\theta_j)
\end{equation}
The observation $y_i$ is distributed according to a mixture, where $\phi = (\phi_1, \ldots, \phi_k)$ are the mixture weights, $\theta = (\theta_1, \ldots, \theta_k)$ are the mixture component parameters, and $p(y_i|\theta_j)$ is the density for the $j$th mixture component. In terms of the original post, $\theta_j = (\alpha_j,\beta_j)$, $p(y_i|\theta_j) = \textsf{Beta}(y_i|\alpha_j,\beta_j)$, and $k=2$.
Assuming independence among the observations conditional on the parameters, the joint likelihood is given by
\begin{equation}
p(y_{1:n}|\phi,\theta) = \prod_{i=1}^n p(y_i|\phi,\theta) ,
\end{equation}
where $y_{1:n} = (y_1, \ldots, y_n)$. The original post suggests a prior for $(\phi,\theta)$ is available. Then Bayes rule delivers
\begin{equation}
p(\phi,\theta|y_{1:n}) \propto p(\phi,\theta)\,p(y_{1:n}|\phi,\theta) . 
\end{equation}
In one sense we are done. In another sense we are just getting going.
An observation by the OP provides a clue as to how to structure the problem so as to make inference more manageable. Let $z_i$ denote a latent "classification" variable, where $z_i \in \{1, \ldots, k\}$. In other words, $z_i = j$ means that $y_i$ comes from component $j$. Of course we don't know for sure the value of $z_i$, but conditional on the parameters the probability that $z_i$ comes from component $j$ is
\begin{equation}
p(z_i=j|\phi,\theta) = \frac{\phi_j\,p(y_i|\theta_j)}{\sum_{j'=1}^k \phi_{j'}\,p(y_i|\theta_{j'})} .
\end{equation}
We can make draws for each classification variable using these probabilities.
Next, conditional on the classifications, we "know" which observations are associated with each component. Let $Y_j = \{y_i|z_i=j\}$. Then for each $j$
\begin{equation}
p(\theta_j|Y_j) \propto p(\theta_j)\,\prod_{y_i\in Y_j} p(y_i|\theta_j) .
\end{equation}
Finally, we can update $\phi$. Let $c = (c_1, \ldots, c_k)$ where $c_j = |Y_j|$ is the number of observations classified with component $j$. Note that $c_j \in \{0, \ldots, n\}$ and $\sum_{j=1}^k c_j = n$. The likelihood for $\phi$ is given by
\begin{equation}
p(c|\phi) = \textsf{Multinomial}(c|n,\phi) .
\end{equation}
If $k = 2$, then $p(c_1|\phi_1) = \textsf{Binomial}(c_1|n,\phi_1)$.
For convenience, assume prior independence between $\phi$ and $\theta$: $p(\phi,\theta) = p(\phi)\,p(\theta)$. Note that $\sum_{j=1}^k \phi_k = 1$ and $\phi_j \ge 0$. As shorthand notation, we can write $\phi \in \Delta^{k-1}$ where $\Delta^{k-1}$ is the $(k-1)$-dimensional simplex. Again for convenience, let us assume
\begin{equation}
p(\phi) = \textsf{Dirichlet}(\phi|\lambda), 
\end{equation}
where $\lambda = (\lambda_1, \ldots, \lambda_k)$. When $k = 2$, the Dirichlet distribution can be expressed as a beta distribution for $\phi_1$. The Dirichlet distribution is the conjugate prior for the multinomial distribution, so the conditional posterior for $\phi$ is also Dirichlet:
\begin{equation}
p(\phi|c) = \textsf{Dirichlet}(\phi|\lambda+c) .
\end{equation}
We now have all the pieces for a Gibbs sampler and we are now in business.
