1 Example of model based clustering with two groups
We will use a toy example to explain the EM algorithm for model based clustering.
Two groups \(k=1,2\)
Univariate observations are observed \(x_i\) with \(i= 1, \dots, n\)
The data of group j follows a normal distribution \(N(u_j,1)\) with a variance \(\sigma^2=1\) and a mean \(\mu_j\) that depends on the group \(j=1,2\).
The data follows the following mixture distribution:
\[f(x) = \pi_1 f_1(x) + (1-\pi_1) f_2(x)\]
With
\(\pi_1\) the probability that a random sample from the population belongs to group 1
\(f_1(x)\) the density of the data in group 1, i.e. \(N(\mu_1,1)\)
\(f_2(x)\) the density of the data in group 2, i.e. \(N(\mu_2,1)\).
Unknowns? Group membership \(z_i\), the group means \(\mu_j\) and the proportion of subjects in the population of group 1 \(\pi_1\) are unknown.
We can estimate the model parameters \(\boldsymbol{\theta} = [\mu_1,\mu_2,\pi_1]^T\) using maximum likelihood.
data.frame(x,zReal) |>ggplot(aes(zReal,x,color=zReal)) +geom_point() +scale_colour_gradient2(low ="blue",mid="white",high="red",midpoint =0.5) +geom_point(x=.95,y=mu1Real, shape=25, col="red", size=3) +geom_point(x=0.05,y=mu2Real, shape=25, col="blue", size=3) +ggtitle("Simulated data with real group means (triangles)")
3 Parameter estimation
If we would know the cluster membership \(z_{i1}\)\[z_{i1}=\left\{ \begin{matrix} 1 & \text{if } x_i \text{ belongs to group 1}\\
0 & \text{if } x_i \text{ belongs to group 2}
\end{matrix}\right.\]
and the joint distribution of \(z_{i1}\) and \(x_i\) then becomes \[\begin{eqnarray}
f(x_i, z_{i1})&=&f(x_i\vert z_{i1}) f(z_{i1})\\
&=& f_1(x_i)^{z_{i1}}f_2(x_i)^{1-z_{i1}}\pi_1^{z_{i1}}(1-\pi_1)^{(1-z_{i1})}
\end{eqnarray}\]
and the log likelihood of the complete data becomes
z <-as.double(x >0)n <-length(z)p <-data.frame(x,zReal,z) |>ggplot(aes(zReal,x,color=z)) +geom_point() +scale_colour_gradient2(low ="blue",mid="white",high="red",midpoint =0.5) +geom_point(x=.95,y=mu1Real, shape=25, col="red", size=3) +geom_point(x=0.05,y=mu2Real, shape=25, col="blue", size=3) +ggtitle("Simulated data with real group means (triangles)")p
4.2 EM algorithm
for (k in1:10){## M-step pi1 <-sum(z)/n mu1 <-sum(z * x)/sum(z) mu2 <-sum((1-z) * x)/sum(1-z)## E-step d1 <-dnorm(x,mean=mu1) d2 <-dnorm(x,mean=mu2) d <- pi1 * d1 + (1-pi1)*d2 z <- pi1*d1/d p <-data.frame(x,zReal,z) |>ggplot(aes(zReal,x,color=z)) +geom_point() +scale_colour_gradient2(low ="blue",mid="white",high="red",midpoint =0.5) +geom_point(x=.95,y=mu1Real, shape=25, col="red", size=3) +geom_point(x=0.05,y=mu2Real, shape=25, col="blue", size=3) +geom_point(x=.95,y=mu1, shape=3, col="red", size=3) +geom_point(x=.05,y=mu2, shape=3, col="blue", size=3) +ggtitle(paste0("iteration ",k, " real mean (triangle), estimate mean (+)"))print(p)}
