EM算法
本文分享了EM算法的两个实践案例，一个是关于两分类的GMM模型，一个是多项分布的混合主题模型。
1.Two-component Gaussian mixture model
1.1
We can know that
\[f(x|\theta)=\tau N(x|\mu_1,\sigma_1^2)+(1-\tau)N(x|\mu_2,\sigma_2^2)\]
\[\begin{split} L(X, Z|\theta) &= \prod_{i=1}^n \pi_k^{Z_{i1}} N(x_i|\theta)^{Z_{i1}}\pi_k^{Z_{i2}} N(x_i|\theta)^{Z_{i2}} \end{split}\]
\[l(\theta|x,z)=\sum_{i=1}^n Z_{i,1} \left( \log (\pi_1) + \log (N(x_i|\mu_1, \sigma_1^2) )\right)+Z_{i,2} \left( \log (\pi_2) + \log (N(x_i|\mu_2, \sigma_2^2) )\right)\]
\[\begin{split} Q(\theta|\theta^{(t)})&=E[\sum_{i=1}^n[Z_{i,1}logN(x_i|\mu_1,\sigma_1^2)+Z_{i,2}logN(x_i|\mu_2,\sigma_2^2)]]+constant\\&:=\sum_{i=1}^n[E[Z_{i,1}|x,\theta^{(t)}][-\frac{1}{2}log(2\pi \sigma^2_1)-\frac{1}{2\sigma_1^2}(x_i-\mu_1)^2]+E[Z_{i,2}|x,\theta^{(t)}][-\frac{1}{2}log(2\pi \sigma^2_2)-\frac{1}{2\sigma_2^2}(x_i-\mu_2)^2]] \end{split}\]
1.2
Do the M-step
\[Max\space \space Q(\theta|\theta^{(t)})\]
\[[\mu_k]:-0.5\sum_{i=1}^nZ_{i,k}^{(t)}(-2)(x_i-\mu_k)=0\] So
\[\mu_k^{(t+1)}=\frac{\sum_{i=1}^nZ_{i,k}^{(t)}x_i}{\sum_{i=1}^nZ_{i,k}^{(t)}}\] \[[\sigma^2_k]:\sum_{i=1}^nZ_{i,k}^{(t)}(-\frac{1}{2 \sigma_k^2}+\frac{(x_i-\mu_k)^2}{2\sigma_k^4})=0\] \[\sigma_k^2{(t+1)}=\frac{\sum_{i=1}^nZ_{i,k}^{(t)}(x_i-\mu_k)^2}{\sum_{i=1}^nZ_{i,k}^{(t)}}\] \[Z_{i,k}^{(t+1)}=\frac{\sum_{i=1}^nZ_{i,k}^{(t)}}{n}\]
1.3
The function called my_gmm_em is the function of GMM model with E-M algorithm.
my_gmm_em=function(k,samp,iteration){
n=length(samp)
prob <- matrix(rep(0, k*n), nrow = n)
weight <- matrix(rep(0, k*n), nrow = n)
threshold <- 1e-5
tau <- c(0.5, 0.5)
miu <- runif(k)
sigma <- runif(k)
out_tau=matrix(rep(0,k*iteration),nrow=iteration)
out_miu=matrix(rep(0,k*iteration),nrow=iteration)
out_sigma=matrix(rep(0,k*iteration),nrow=iteration)
out_tau[1,]=tau
out_miu[1,]=miu
out_sigma[1,]=sigma
for (step in 1:iteration) {
for (j in 1:k) {
prob[, j] <- sapply(samp, dnorm, miu[j], sigma[j])
weight[, j] <- tau[j] * prob[, j]
}
row_sum <- rowSums(weight)
prob <- weight/row_sum
oldtau <- tau
oldmiu <- miu
oldsigma <- sigma
for (j in 1:k) {
sum1 <- sum(prob[, j])
sum2 <- sum(samp*prob[, j])
tau[j] <- sum1/n
miu[j] <- sum2/sum1
sum3 <- sum(prob[, j]*(samp-miu[j])^2)
sigma[j] <- sqrt(sum3/sum1)
}
out_tau[step,]=tau
out_miu[step,]=miu
out_sigma[step,]=sigma
#if (sum(abs(tau - oldtau)) < threshold &
#sum(abs(miu - oldmiu)) < threshold &
#sum(abs(sigma - oldsigma)) < threshold) break
}
return(mylist=list("tau"=out_tau,"miu"=out_miu,"sigma"=out_sigma))
}The code below aims to generate the data, and then run the GMM code. The out put is \(\theta\)
set.seed(100000)
n <- 100
tau1 <- 0.3
miu1 <- 0
sigma1 <- 1
tau2 <- 0.7
miu2 <- 2
sigma2 <- 0.5
n1 <- floor(n*tau1)
n2 <- n-n1
samp <-numeric(n)
samp[1:n1] <- rnorm(n1, miu1, sigma1)
samp[(n1+1):n] <- rnorm(n2, miu2, sigma2)
out=my_gmm_em(2,samp,200)print(paste("Final Estimate of tau = (",out$tau[200,1],",",out$tau[200,2],")"))## [1] "Final Estimate of tau = ( 0.323950241780784 , 0.676049758219216 )"print(paste("Final Estimate of miu = (",out$miu[200,1],",",out$miu[200,2],")"))## [1] "Final Estimate of miu = ( 0.0606505669653436 , 1.99868213429268 )"print(paste("Final Estimate of sigma = (",out$sigma[200,1],",",out$sigma[200,2],")"))## [1] "Final Estimate of sigma = ( 1.03062008863022 , 0.473597160166307 )"The output above shows the estimation of \(\tau\), \(\mu\), \(\sigma\). They are similar to the real value of \(\tau\), \(\mu\), \(\sigma\).
As we can find that, the size of the data set is just 200, it’s a little bit small for GMM model, the larger the size is, more accurate the estimation will be.
1.4
When it comes to the convergence of the EM algorithm, we can plot the list of the output estimation.
I calculate the error of the \(\tau\), \(\mu\), \(\sigma\) in the way shown below, take \(\tau\) as an example
\[err_\tau=\sum_{i=1}^k(\hat{\tau_i}-\tau_i)^2\]
err_tau=rowSums((sweep(out$tau,2,c(0.3,0.7)))^2)
plot(err_tau, main="Error of the Estimation of tau")
err_miu=rowSums((sweep(out$miu,2,c(0,2)))^2)
plot(err_miu, main="Error of the Estimation of mu")
err_sigma=rowSums((sweep(out$sigma,2,c(1,0.5)))^2)
plot(err_sigma, main="Error of the Estimation of sigma")
From the picture above, we can see the convergence of the six parameters.
2.A (three-category) multinomial mixture problem
2.1 EM steps
The multinomial random variable with the following probability mass function
\[f(x|\theta)=\frac{(x_{1}+x_2+x_3)!}{x_{1}!x_2!x_3!}(\frac{4+\theta}{6})^{x_{1}}(\frac{1-\theta}{3})^{x_2}(\frac{\theta}{6})^{x_3}\]
We have
\(p_1=\frac{4+\theta}{6}\), \(p_2=\frac{1-\theta}{3}\), \(p_3=\frac{\theta}{6}\)
We know that
\[\sum_{i=1}^3p_i=1\] we can split \(p_1\) into two parts \(p_{11}=\frac{2}{3}\), \(p_{12}=\frac{\theta}{6}\), so we can get the hidden variable \(p_{12}\)
\[f(x|\theta)=\frac{(x_{11}+x_{12}+x_2+x_3)!}{x_{11}!x_{12!}!x_2!x_3!}(\frac{2}{3})^{x_{11}}(\frac{\theta}{6})^{x_{12}}(\frac{1-\theta}{3})^{x_2}(\frac{\theta}{6})^{x_3}\]
\[l(\theta|X)=ln(\frac{(x_{11}+x_{12}+x_2+x_3)!}{x_{11}!x_{12!}!x_2!x_3!})+x_{11}ln(\frac{2}{3})+x_{12}ln(\frac{\theta}{6})+x_{2}ln(\frac{1-\theta}{3})+x_{3}ln(\frac{\theta}{6})\]
\[Q(\theta|\theta^{(t+1)})=E_{\theta}[ln(\frac{(x_{11}+x_{12}+x_2+x_3)!}{x_{11}!x_{12!}!x_2!x_3!})+x_{11}ln(\frac{2}{3})+x_{12}ln(\frac{\theta}{6})+x_{2}ln(\frac{1-\theta}{3})+x_{3}ln(\frac{\theta}{6})]\]
The expectation is about \(\theta\), so we can have
\[Q(\theta|\theta^{(t+1)})=ln(\frac{(x_{11}+x_{12}+x_2+x_3)!}{x_{11}!x_{12!}!x_2!x_3!})+E_{\theta}(x_{11})ln(\frac{2}{3})+E_{\theta}(x_{12})ln(\frac{\theta}{6})+E_{\theta}(x_{2})ln(\frac{1-\theta}{3})+E_{\theta}(x_{3})ln(\frac{\theta}{6})\]
So when we remove the constant value
\[Q(\theta|\theta^{(t+1)})=E_{\theta}(x_{12}|x_1)ln(\frac{\theta}{6})+x_{2}ln(\frac{1-\theta}{3})+x_{3}ln(\frac{\theta}{6})\]
We know that
\[E_{\theta}(x_{12}|x_1)=\frac{\frac{\theta}{6}}{\frac{4+\theta}{6}}x_1=\frac{\theta}{4+\theta}x_1\]
we can denote it as \[E_{\theta^{(0)}}(x_{12}|x_1)=\frac{\theta^{(0)}}{4+\theta^{(0)}}x_1=x_{12}^{(0)}\] So we have Q function that
\[Q(\theta|\theta^{(t+1)})=x_{12}^{(t+1)}ln(\frac{\theta}{6})+x_{2}ln(\frac{1-\theta^{(t+1)}}{3})+x_{3}ln(\frac{\theta^{(t+1)}}{6})\]
Above is the E-step. To maximum Q function above, we can have M-step below
\[[\theta^{(t+1)}]:\frac{6 x_{12}^{(t+1)}}{\theta^{(t+1)}}-\frac{3x_{2}}{1-\theta^{(t+1)}}+\frac{6x_{3}}{\theta^{(t+1)}}=0\] We finally have
\[\theta^{(t+1)}=\frac{2x_{12}^{(t)}+2x_3}{2x_{12}^{(t)}+x_2+2x_3}\]
\[x_{12}^{(t)}=\frac{\theta^{(t)}}{4+\theta^{(t)}}x_1\]
2.2
my_multno_em=function(samp,iteration){
n=length(samp)
threshold <- 1e-5
theta_ini = 0.5
x12_ini=(theta_ini/(4+theta_ini))*samp[1]
out_theta=matrix(rep(0,iteration),nrow=iteration)
out_x12=matrix(rep(0,iteration),nrow=iteration)
out_theta[1,]=theta_ini
out_x12[1,]=x12_ini
for (step in 2:iteration) {
out_x12[step]=(out_theta[step-1]/(4+out_theta[step-1]))*samp[1]
out_theta[step]=(2*out_x12[step-1]+2*samp[3])/(2*out_x12[step-1]+2*samp[3]+samp[2])
}
return(mylist=list("theta"=out_theta,"x_12"=out_x12))
}samp=c(42,10,15)
out=my_multno_em(samp,20)plot(out$theta,main="The estimation of theta")
out$theta[20]## [1] 0.8155887Above we can find the estimation of \(\theta\), It converges to 0.8155887.