# 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?

• Take a look at the EM algorithm. For example, Bishop's book (Pattern Recognition and Machine Learning) has a nice chapter on mixture models. – jpmuc Jun 4 '17 at 16:20
• @jpmuc thanks for the comment, I'll look up the reference tomorrow. In the meantime, does the EM algorithm allow for updating information on a "prior" model or is it only used for estimating parameters of a mixture model "from scratch"? – Constantinos Jun 4 '17 at 18:37
• EM is a good suggestion, but if the data are truly streaming (observations arriving sequentially), then you will need an online flavor of EM. I'd suggest looking at arxiv.org/abs/0712.4273. Variational Bayesian inference is also an option. – Nate Pope Jun 4 '17 at 18:40
• @NatePope thanks for that. Would it make any difference if the $p_i$ weren't sequentially arriving but in a batch of $n$ data points? – Constantinos Jun 4 '17 at 19:02
• Nope... If the point is that you want to update current estimates based on incoming data without recomputing everything from scratch, it's still an online learning problem, and the method I linked to is still applicable. – Nate Pope Jun 4 '17 at 19:17

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: $$$$p(y_i|\phi,\theta) = \sum_{j=1}^k \phi_j\,p(y_i|\theta_j)$$$$ 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 $$$$p(y_{1:n}|\phi,\theta) = \prod_{i=1}^n p(y_i|\phi,\theta) ,$$$$ where $$y_{1:n} = (y_1, \ldots, y_n)$$. The original post suggests a prior for $$(\phi,\theta)$$ is available. Then Bayes rule delivers $$$$p(\phi,\theta|y_{1:n}) \propto p(\phi,\theta)\,p(y_{1:n}|\phi,\theta) .$$$$ 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 $$$$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'})} .$$$$ 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$$ $$$$p(\theta_j|Y_j) \propto p(\theta_j)\,\prod_{y_i\in Y_j} p(y_i|\theta_j) .$$$$

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 $$$$p(c|\phi) = \textsf{Multinomial}(c|n,\phi) .$$$$ 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 $$$$p(\phi) = \textsf{Dirichlet}(\phi|\lambda),$$$$ 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: $$$$p(\phi|c) = \textsf{Dirichlet}(\phi|\lambda+c) .$$$$

We now have all the pieces for a Gibbs sampler and we are now in business.