\chapter{Algoritmo EM}

\section{Estimaci\'{o}n de m\'{a}xima verosimilitud.}

Deseamos encontrar para un conjunto de observaciones $x_{i},$\ los
parámetros de la distribución que las genero. Considerando el caso
de una distribución normal, los parámetros $\theta $ a calcular
serán la media $\mu $ y la varianza $\sigma ^{2}$
\begin{equation}
g(x_{i},\mu ,\sigma )=\frac{1}{\sqrt{2\pi }\sigma }\exp \left( \frac{%
-(x_{i}-\mu )^{2}}{2\sigma ^{2}}\right)
\end{equation}

La probabilidad conjunta de las observaciones, asumiendo
independencia, dados los parámetros de la distribución se pueden
calcular como
\[
P(x|\mu ,\sigma )=\prod\limits_{i=1}^{N}g(x_{i},\mu ,\sigma )
\]
y aplicando regla de Bayes podemos calcular la verosimilitud como
\[
P(\mu ,\sigma |x)=\frac{P(x|\mu ,\sigma )P(\mu )P(\sigma )}{P(x)}=P(x|\mu
,\sigma )
\]
considerando $P(\mu ),$ $P(\sigma )$ y $P(x)$ constantes

La logverosimilitud para las $N$ observaciones es:
\begin{eqnarray}
P(\mu ,\sigma |x) &=&\prod\limits_{i=1}^{N}g(x_{i},\mu ,\sigma )  \nonumber
\\
£(\mu ,\sigma ) &=&\ln P(\mu ,\sigma |x)=\sum\limits_{i=1}^{N}\left[ -\ln
\sqrt{2\pi }-\ln \sigma -\frac{(x_{i}-\mu )^{2}}{2\sigma ^{2}}\right]
\end{eqnarray}
El estimador $\theta ^{*}$ de máxima verosimilitud lo obtenemos calculando $%
\frac{\partial £(\theta )}{\partial \theta }=0$%
\begin{eqnarray}
\frac{\partial £(\mu ,\sigma )}{\partial \mu } &=&\sum\limits_{i=1}^{N}%
\left[ \frac{x_{i}-\mu }{\sigma ^{2}}\right] =0  \nonumber \\
\mu ^{*} &=&\frac{\sum\limits_{i=1}^{N}x_{i}}{N}
\end{eqnarray}
\begin{eqnarray}
\frac{\partial £(\mu ,\sigma )}{\partial \sigma } &=&\sum\limits_{i=1}^{N}%
\left[ -\frac{1}{\sigma }+\frac{(x_{i}-\mu )^{2}}{2\sigma ^{3}}\right] =0
\nonumber \\
\left( \sigma ^{2}\right) ^{*} &=&\frac{\sum\limits_{i=1}^{N}(x_{i}-\mu )^{2}%
}{N}
\end{eqnarray}

\section{El Algoritmo EM.}

El t\'{e}rmino EM fue introducido por Dempster en 1977 para designar un
algoritmo muy general para la estimaci\'{o}n de m\'{a}xima verosimilitud en
conjuntos de datos incompletos. El nombre se refiere a los pasos que
constituyen el algoritmo: un paso E de c\'{a}lculo de la esperanza, seguido
por otro paso M, de maximizaci\'{o}n de la verosimilitud. Los pasos a seguir
para el algoritmo EM son:

1.- Se asigna un valor inicial a los par\'{a}metros por alg\'{u}n m\'{e}todo
razonable.

2.- Los valores que faltan son reemplazados por sus valores
esperados, condicionados a los valores observados y a los valores
corrientes de los par\'{a}metros.

3.- Los par\'{a}metros son reestimados usando los valores observados y
modificando los valores faltantes.

4.- Los valores perdidos son recalculados asumiendo que los nuevos
par\'{a}metros son correctos.

5.- Los par\'{a}metros son recalculados una vez m\'{a}s en el paso tres y en
el paso cuatro hasta cumplir cierto criterio de paro.

\section{Aplicación del EM a mezclas de distribuciones.}


El problema cuando tratamos con mezclas de distribuciones es que
uno no sabe que distribución en la mezcla fue la que produjo cada
ejemplo de observación. Si tuviéramos esa información disponible,
solo resta calcular el estimador de máxima verosimilitud para una
distribución en particular. Utilizaremos una variable $\varepsilon
$, para etiquetar la i-esima observación generada por la k-esima
distribución así:
\begin{eqnarray}
\varepsilon _{ik} &=&1\text{ si }x_{i}\thicksim g(x_{i},\theta _{k})
\nonumber \\
\text{ }\varepsilon _{ik} &=&0\text{ para los demás}
\end{eqnarray}
Nuestros datos completos son $X=(x_{i},\varepsilon _{i})$ donde
$x_{i}$ esta dado y$\ \varepsilon _{i}$\ lo desconocemos. La
verosimilitud para los datos completos es $P(\varepsilon ,\theta
|x)=P(x|\varepsilon ,\theta )*P(\varepsilon |\theta )*$ $P(\theta
)/P(x)$\ donde cada término lo definidos como
\begin{equation}
P(x|\varepsilon ,\theta )=\prod\limits_{i=1}^{N}g(x_{i},\theta
_{k})=\prod\limits_{i=1}^{N}\prod\limits_{k=1}^{C}g(x_{i},\theta
_{k})^{\varepsilon _{ik}}
\end{equation}
\begin{equation}
P(\varepsilon |\theta )=\prod\limits_{k=1}^{C}\alpha
_{k}=\prod\limits_{i=1}^{N}\prod\limits_{k=1}^{C}\alpha
_{k}^{\varepsilon _{ik}}
\end{equation}
y $P(\theta )$, $P(x)$ como constantes.

La logverosimilitud $£(\varepsilon ,\theta )=\ln P(\varepsilon ,\theta |x)$
queda definida
\begin{eqnarray}
£(\varepsilon ,\theta ) &=&\ln
\prod\limits_{i=1}^{N}\prod\limits_{k=1}^{C}\left[ \alpha _{k}g(x_{i},\theta
_{k})\right] ^{\varepsilon _{ik}}  \nonumber \\
&=&\sum\limits_{i=1}^{N}\sum\limits_{k=1}^{C}\varepsilon _{ik}\left[ \ln
\alpha _{k}+\ln g(x_{i},\theta _{k})\right]   \label{log}
\end{eqnarray}
Tendremos que determinar dadas las observaciones $x_{i}$ los parámetros $%
\theta _{k}$ de las $C$ distribuciones que conforman la mezcla, así como el
valor de $\alpha _{k}$\ que indica la proporción de observaciones que
pertenecen a la $k-esima$ distribución $g(x_{i,}\theta _{k})$. Como
desconocemos el valor $\varepsilon _{ik}$ utilizaremos el valor esperado $%
\overline{\varepsilon }_{ik}=E\left[ \varepsilon _{ik}|x_{i},\theta \right] $
(Paso E \ del algoritmo EM). Analizando los valores de la variable $%
\varepsilon _{ik}$ ($\varepsilon _{i}$ es un vector de $C-1$ de
ceros y únicamente un $1$) por lo que el valor esperado
$\overline{\varepsilon }_{ik} $ es la probabilidad marginal a
posteriori de que $\varepsilon _{ik}=1$ dadas $x_{i}$ y $\theta $.
Utilizando la regla de Bayes podemos encontrar que:
\begin{eqnarray}
\overline{\varepsilon }_{ik}^{t} &=&P(\varepsilon _{ik}=1|x_{i},\theta )^{t}=%
\frac{P(x_{i}|\varepsilon _{ik}=1,\theta )^{t}P(\varepsilon _{ik}=1|\theta
)^{t}}{P(x_{i}|\theta )^{t}}  \nonumber \\
\overline{\varepsilon }_{ik}^{t} &=&P(\varepsilon _{ik}=1|x_{i},\theta )^{t}=%
\frac{g(x_{i},\theta _{k}^{t})*\alpha _{k}^{t}}{\sum%
\limits_{k=1}^{C}g(x_{i},\theta _{k}^{t})*\alpha _{k}^{t}}  \label{iii}
\end{eqnarray}
Para calcular (\ref{iii}) tenemos que maximizar la log-verosimilitud de las
observaciones sujeta a la restricción de que $\sum\limits_{k}\alpha _{k}=1$.
Aplicando multiplicadores de Lagrange tenemos:
\begin{equation}
l_{obs}=\sum\limits_{i=1}^{N}\sum\limits_{k=1}^{C}\overline{\varepsilon }%
_{ik}\left[ \ln \alpha _{k}+\ln g(x_{i},\theta _{k})\right] -\lambda \left[
\sum\limits_{k=1}^{C}\ p_{k}-1\right]   \label{iv}
\end{equation}
Para resolver calculamos $\frac{\partial l_{obs}}{\partial \alpha _{k}}=0$
para la ecuación (\ref{iv})
\begin{equation}
\frac{\partial l_{obs}}{\partial \alpha _{k}}=\sum\limits_{i=1}^{N}\frac{%
\overline{\varepsilon }_{ik}}{\alpha _{k}}-\lambda =0  \label{v}
\end{equation}
Reorganizando términos y sumando sobre las $C$\ distribuciones tendremos
\begin{eqnarray}
\lambda \sum\limits_{k=1}^{C}\alpha _{k}
&=&\sum\limits_{i=1}^{N}\sum\limits_{k=1}^{C}\overline{\varepsilon }_{ik}
\nonumber \\
\lambda  &=&N  \label{vi}
\end{eqnarray}
Sustituyendo en la ecuación (\ref{v}) el valor de $\lambda $%
\begin{equation}
\alpha _{k}=\frac{\sum\limits_{i=1}^{N}\overline{\varepsilon }_{ik}}{N}
\end{equation}

\subsection{Mezclas finitas de distribuciones normales.}

Para el desarrollo del algoritmo, ser\'{a} necesario conocer a priori el
tipo distribuciones que forman la mezcla. Un caso sencillo resulta ser
cuando utilizamos una distribuci\'{o}n normal. El vector de parámetros $%
\theta _{k}$ es $(\mu _{k},\sigma _{k})$ y la ditribuci\'{o}n esta dada por:
\begin{equation}
g(x_{i},\mu _{k},\sigma _{k})=\frac{1}{\sqrt{2\pi }\sigma _{k}}\exp \left(
\frac{-(x_{i}-\mu _{k})^{2}}{2\sigma _{k}^{2}}\right)
\end{equation}
Dada la funci\'{o}n de distribuci\'{o}n la log-verosimilitud es:
\begin{equation}
\pounds (\mu ,\sigma )=\sum\limits_{i=1}^{N}\sum\limits_{k=1}^{C}\overline{%
\varepsilon }_{ik}\left[ \ln \alpha _{k}-\ln \sqrt{2\pi }-\ln \sigma -\frac{%
(x_{i}-\mu )^{2}}{2\sigma ^{2}}\right]
\end{equation}
Maximizamos $\pounds (\mu ,\sigma )$ con respecto al vector de
par\'{a}metros $(\mu ,\sigma )$%
\begin{equation}
\frac{\partial E}{\partial \mu _{k}}=\sum\limits_{i=1}^{N}\overline{%
\varepsilon }_{ik}\frac{(x_{i}-\mu _{k})}{\sigma _{k}^{2}}=0
\end{equation}
Resolviendo para $\mu _{k}$ y recordando que es el valor a tiempo $t+1.$%
\begin{equation}
\left( \mu _{k}^{\ast }\right) ^{t+1}=\frac{\sum\limits_{i=1}^{N}\overline{%
\varepsilon }_{ik}^{t}x_{i}}{\sum\limits_{i=1}^{N}\overline{\varepsilon }%
_{ik}^{t}}
\end{equation}
Hacemos lo mismo para el caso de $\sigma _{k}$ .
\begin{equation}
\frac{\partial E}{\partial \sigma _{k}}=\sum\limits_{i=1}^{N}\overline{%
\varepsilon }_{ik}\left( -\frac{1}{\sigma _{k}}+\frac{(x_{i}-\mu _{k})^{2}}{%
\sigma _{k}^{3}}\right) =0
\end{equation}

\begin{equation}
\left( \sigma _{k}^{\ast 2}\right) ^{t+1}=\frac{\sum\limits_{i=1}^{N}%
\overline{\varepsilon }_{ik}^{t}(x_{i}-\mu _{k})^{2}}{\sum\limits_{i=1}^{N}%
\overline{\varepsilon }_{ik}^{t}}
\end{equation}

\section{Imágenes con tonos de gris constantes a pedazos.}

Para este tipo de imágenes queremos determinar el valor del tono
de gris que caracteriza a una region $R_{k}$ tal que $I(r)=\mu
_{k}+n_{r}$ $\forall r\in R_{k}.$ A priori sabemos que los tonos
de gris son constantes a pedazos y que existe un ruido de
observaci\'{o}n $n_{r}$ que es normal $N(0,\sigma ^{2})$. La
verosimilitud de los datos completos esta dada por

\begin{equation}
P(x,\varepsilon |\theta )=P(x|\varepsilon ,\theta )\ast P(\varepsilon
|\theta )
\end{equation}
El t\'{e}rmino $P(x|\varepsilon ,\theta )$ es la probabilidad de una
observaci\'{o}n $x$ dados $\varepsilon ,\theta $. Consideraremos que la
distribuci\'{o}n es gaussiana y que el tono de gris a calcular tiene media $%
\mu _{k}$ y varianza $\sigma _{k}^{2}=1$.

\begin{equation}
P(x|\varepsilon ,\theta
)=\prod\limits_{i=1}^{N}\prod\limits_{k=1}^{C}g(x_{i},\mu
_{k},1)^{\varepsilon _{ik}}
\end{equation}
El modelo m\'{a}s popular, usado para restringir un campo a ser
constante a pedazos, es el modelo generalizado de Ising, cuyos
potenciales están dados por.
\begin{eqnarray}
V(\varepsilon _{i},\varepsilon _{j}) &=&-1\ si\ \varepsilon _{i}=\varepsilon
_{k}  \nonumber \\
V(\varepsilon _{i},\varepsilon _{j}) &=&+1\ si\ \varepsilon _{i}\neq
\varepsilon _{j}
\end{eqnarray}
La probabilidad conjunta $P(\varepsilon |\theta )$ no prodr\'{a}
ser representada como el producto de las probabilidades
marginales, ya que el imponer una restricci\'{o}n de coherencia
espacial hace que se pierda la independencia. Supondremos una
distribuci\'{o}n de probabilidades de Gibbs para $P(\varepsilon
|\theta )$ (campo aleatorio Markoviano) dada por.
\begin{equation}
P(\varepsilon |\theta )=e^{-\beta \sum\limits_{<i,j>}V(\varepsilon
_{i},\varepsilon _{j})}
\end{equation}
Las expresiones finales para la verosimilitud y la logverosimilutud son:
\begin{equation}
P(x,\varepsilon |\theta
)=\prod\limits_{i=1}^{N}\prod\limits_{k=1}^{C}g(x_{i},\mu
_{k},1)^{\varepsilon _{ik}}\times e^{-\beta \sum\limits_{<i,j>}V(\varepsilon
_{ik},\varepsilon _{jk})}
\end{equation}
\begin{equation}
\pounds (\varepsilon ,\mu )=-\sum\limits_{i=1}^{N}\sum\limits_{k=1}^{C}\frac{%
(x_{i}-\mu _{k})^{2}}{2}\varepsilon _{ik}-\beta
\sum\limits_{<i,j>}V(\varepsilon _{ik},\varepsilon _{jk})
\end{equation}
Minimizar la expresión anterior resulta complicada, para ello se
utilizan algoritmos como Simulated Annealing, el muestreador
Gibbsiano, el algoritmo de Montecarlo, y de manera mas eficaz
campos de medidas.
