Algoritmo EM para MLE de una muestra normal bivariada con datos perdidos: Atascado en el paso M

Question

Estoy tratando de entender la aplicación del algoritmo EM para calcular el MLE en un problema de datos perdidos.

En concreto, supongamos que $(x_1,y_1),\ldots,(x_n,y_n)$ es una muestra aleatoria de la distribución normal bivariada con media $(0,0)$ y la matriz de covarianza desconocida $$ \Sigma=\begin{bmatrix} \sigma_x^2&\rho\sigma_x\sigma_y\\ \rho\sigma_x\sigma_y&\sigma_y^2 \end{bmatrix}. $$ Quiero encontrar la MLE de $\Sigma$ dado que el primer $\nu< n$ de la $y$ -falta de coordenadas.

Para el paso E del algoritmo EM, necesito calcular

$$ Q_i(\Sigma,\tilde\Sigma):= \begin{cases} \mathbb{E}_{y\sim f(y\mid x_i,\tilde\Sigma)}[\log f(x_i,y\mid\Sigma)]& \text{if }1\leq i\leq \nu,\\ \log f(x_i,y_i\mid\Sigma)&\text{if }\nu< i\leq n, \end{cases} $$ donde $$ \begin{aligned} -\log f(x,y\mid\Sigma) = \log2\pi &+ \frac12\log\sigma_x^2 + \frac12\log\sigma_y^2 + \frac12\log(1-\rho^2)\\ &+ \frac{1}{2(1-\rho^2)} \left\{\frac{x^2}{\sigma_x^2} - \frac{2\rho xy}{\sigma_x\sigma_y} + \frac{y^2}{\sigma_y^2}\right\} \end{aligned} $$ y $$ \tilde\Sigma=\begin{bmatrix} \tilde\sigma_x^2&\tilde\rho\tilde\sigma_x\tilde\sigma_y\\ \tilde\rho\tilde\sigma_x\tilde\sigma_y&\tilde\sigma_y^2 \end{bmatrix}. $$

Desde $$ \begin{aligned} \mathbb{E}_{y\sim f(y\mid x_i,\tilde\Sigma)} [y] &= \frac{\tilde\rho \tilde\sigma_yx_i}{\tilde\sigma_x},\\ \mathbb{E}_{y\sim f(y\mid x_i,\tilde\Sigma)}[y^2] &=\tilde\sigma_y^2(1-\tilde\rho^2) + \frac{\tilde\rho^2 \tilde\sigma^2_yx_i^2}{\tilde\sigma_x^2}, \end{aligned} $$

Lo entiendo.

$$ \begin{aligned} -Q_i(\Sigma,\tilde\Sigma) = \log2\pi &+ \frac12\log\sigma_x^2 + \frac12\log\sigma_y^2 + \frac12\log(1-\rho^2)\\ &+ \frac{1}{2(1-\rho^2)} \left\{\frac{x_i^2}{\sigma_x^2} - \frac{2\rho\tilde\rho\tilde\sigma_y^2 x_i^2}{\sigma_x\tilde\sigma_x\sigma_y} + \frac{\tilde\sigma_x^2\tilde\sigma_y^2(1-\tilde\rho^2)+ \tilde\rho^2 \tilde\sigma^2_yx_i^2}{\tilde\sigma_x^2\sigma_y^2} \right\}. \end{aligned} $$

para $1\leq i\leq \nu$ .

Ahora, para el paso M, necesito calcular $$ \operatorname*{argmax}_{\Sigma} \sum_{i=1}^n Q_i(\Sigma,\tilde\Sigma). $$

Y aquí estoy atascado. ¿Existe una forma agradable para el óptimo $\Sigma$ ?

Answer 1

Una pista: Al considerar \begin{aligned} \log f(x_i,y_i\mid\Sigma) = -\log2\pi &- \frac12\log\sigma_x^2 - \frac12\log\sigma_y^2 - \frac12\log(1-\rho^2)\\ &+ \frac{1}{2(1-\rho^2)} \left\{\frac{x_i^2}{\sigma_x^2} - \frac{2\rho x_iy_i}{\sigma_x\sigma_y} + \frac{y^2}{\sigma_y^2}\right\} \end{aligned} y \begin{aligned} Q_i(\Sigma,\tilde\Sigma) = &-\log2\pi- \frac12\log\sigma_x^2 - \frac12\log\sigma_y^2 - \frac12\log(1-\rho^2)\\ &+ \frac{1/2}{(1-\rho^2)} \left\{\frac{x_i^2}{\sigma_x^2} - \frac{2\rho x_i\tilde\rho\tilde\sigma_y x_i}{\sigma_x\tilde\sigma_x\sigma_y} + \frac{\tilde\sigma_y^2(1-\tilde\rho^2)+\tilde\rho^2\tilde\sigma_y^2 x_i^2/\tilde\sigma_x^2}{\sigma_y^2} \right\} \end{aligned} ambas expresiones son esencialmente de formas idénticas como funciones de $\Sigma$ . Esto significa que la función objetivo a optimizar se escribe como \begin{align}\sum_{i=1}^\nu\, &\log f(x_i,\mathbb E_{y\sim f(y\mid x_i,\tilde\Sigma)}[y_i]\mid\Sigma)+\\ \sum_{i=1}^\nu\, &\left\{\log f(0,\tilde\sigma_y(1-\tilde\rho^2)^{1/2})\mid\Sigma)+\log2\pi + \frac12\log[\sigma_x^2 \sigma_y^2(1-\rho^2)]\right\}+\\\sum_{i\nu+=1}^n\, &\log f(x_i,\tilde y_i\mid\Sigma)\end{align} es decir, como normal log-likelihood para una muestra normal modificada $\mathbf Z$ (dependiendo de la corriente $\tilde\Sigma$ ) $$ -\frac n2 \log|\Sigma|-\frac12\text{trace}(\mathbf Z \Sigma^{-1}) $$ El estimador de $\Sigma$ por lo tanto, se puede derivar como en el caso normal .

Answer 2

Para completar la respuesta de @xian (¡y para asegurarme de que la entiendo!):

Para $1\leq i\leq \nu$ , dejemos que $$ \begin{aligned} \tilde y_i & = \mathbb{E}_{y\sim f(y\mid x_i,\tilde\Sigma)}[y],\\ \tilde z_i^2 & = \operatorname{Var}_{y\sim f(y\mid x_i,\tilde\Sigma)}[y]. \end{aligned} $$

Lo tenemos: $$ \tilde z_i^2 = \mathbb{E}_{y\sim f(y\mid x_i,\tilde\Sigma)}[y^2] - \tilde y_i^2 $$

Si $1\leq i\leq \nu$ , $$ \begin{aligned} Q_i(\Sigma,\tilde\Sigma) &= \mathbb{E}_{y\sim f(y\mid x_i,\tilde\Sigma)}[\log f(x_i,y\mid\Sigma)]\\ &= \mathbb{E}_{y\sim f(y\mid x_i,\tilde\Sigma)}\left[-\log2\pi - \frac12\log\Sigma - \frac{1}{2(1-\rho^2)} \left\{\frac{x_i^2}{\sigma_x^2} - \frac{2\rho x_iy}{\sigma_x\sigma_y} + \frac{y^2}{\sigma_y^2}\right\}\right]\\ &= -\log2\pi - \frac12\log\Sigma - \frac{1}{2(1-\rho^2)} \left\{\frac{x_i^2}{\sigma_x^2} - \frac{2\rho x_i\tilde y_i}{\sigma_x\sigma_y} + \frac{\mathbb{E}_{y\sim f(y\mid x_i,\tilde\Sigma)}[y^2]}{\sigma_y^2}\right\}\\ &= -\log2\pi - \frac12\log\Sigma - \frac{1}{2(1-\rho^2)} \left\{\frac{x_i^2}{\sigma_x^2} - \frac{2\rho x_i\tilde y_i}{\sigma_x\sigma_y} + \frac{\tilde y_i^2}{\sigma_y^2}\right\} - \frac{\tilde z_i^2}{2\sigma_y^2(1-\rho^2)}\\ &= \log f(x_i,\tilde y_i\mid\Sigma) - \frac{\tilde z_i^2}{2\sigma_y^2(1-\rho^2)}. \end{aligned} $$

El término $\log f(x_i,\tilde y_i\mid\Sigma)$ es intuitivamente atractiva; hemos imputado la incógnita $y$ -valores con sus expectativas condicionales.

Consideremos el término de corrección:

$$ - \frac{\tilde z_i}{2\sigma_y^2(1-\rho^2)} = \log f(0,\tilde z_i\mid\Sigma) + \frac12\log 2\pi + \frac12\log\det\Sigma. $$

Así, si $1\leq i\leq \nu$ , $$ Q_i(\Sigma,\tilde\Sigma) = \log f(x_i,\tilde y_i\mid\Sigma) + \log f(0,\tilde z_i\mid\Sigma) + \frac12\log 2\pi + \frac12\log\det\Sigma. $$

El paso M consiste en maximizar la función $$ \begin{aligned} Q(\Sigma,\tilde\Sigma) &:= \sum_{i=1}^n Q_i(\Sigma,\tilde\Sigma) \end{aligned} $$

Al ser esencialmente la suma de las probabilidades logarítmicas gaussianas y los términos de corrección "simples", podemos maximizar $Q(\Sigma,\tilde\Sigma)$ utilizando técnicas estándar. A no ser que haya metido la pata (probablemente), el maximizador resultante es

$$ \Sigma = \frac1n\begin{bmatrix} \sum_{i=1}^n x_i^2& \sum_{i=1}^\nu x_i \tilde y_i + \sum_{i=\nu+1}^n x_i y_i\\ \sum_{i=1}^\nu x_i \tilde y_i + \sum_{i=\nu+1}^n x_i y_i& \sum_{i=1}^\nu (\tilde y_i^2 + \tilde{z}_i^2) + \sum_{i=\nu+1}^n y_i^2 \end{bmatrix}. $$