Estadística computacional

Tema 1: El algoritmo EM
Sesión 4
Doctorado en Estadística

Profesor: Dr. Harvey Rosas

Alumno: Marcelo Rodríguez

Facultad de Ciencias
Universidad de Valparaíso

Valparaíso, Chile, 23 de abril de 2015

Contenidos

1 Introduccion

2 Un ejemplo inspirador

3 El algortitmo EM

4 Ejemplos del algortitmo EM

Introduccion
 I II III IV El algoritmo EM

Introducción

El algoritmo EM (Expectation-Maximization) es una tecnica de

optimizacion originalmente introducida por Dempster, Laird and
Rubin (1977), en su publicación Maximum Likelihood from In-
complete Data via the EM Algorithm. Se utiliza en estadística
para encontrar estimadores de verosimilitud maxima (VM) de pa-
rámetros en modelos probabilísticos que dependen de variables no
observables (datos perdidos).

El algoritmo EM alterna pasos de esperanza (paso E), don-

de se calcula la esperanza de la verosimilitud mediante la
inclusión de variables latentes como si fueran observables,

y un paso de maximización (paso M), donde se calculan esti-

madores de VM de los parámetros mediante la maximización
de la verosimilitud esperada del paso E.

Los parámetros que se encuentran en el paso M se usan para

comenzar el paso E siguiente, y así el proceso se repite.

[email protected] IEUV Estadística computacional 15-10-2014 4/25

Un ejemplo inspirador
 I II III IV El algoritmo EM

Ejemplo inspirador: estimación de ligamiento

en genética (Rao, 1973)

Este problema es propuesto por Rao (1973, pp. 368-369) y ana-

lizado por Dempester (1977) para presentar el algoritmo EM.
Se presenta datos en que 197 animales se distribuyen en cuatro
categorías (AB, Ab, aB y ab), de modo que el vector de datos
observados de frecuencias es

x = (x1 , x2 , x3 , x4 )T = (125, 18, 20, 34)T .

Se postula que provienen de una distribucion multinomial con
cuatro categorias con funcion de probabilidad

Pr(X1 = x1 , X2 = x2 , X3 = x3 , X4 = x4 ) = g(x|θ),
es decir
 x1  x2  x3   x4
n! 1 1 1 1 1 1 1
g(x|θ) = + θ − θ − θ θ .
x1 !x2 !x3 !x4 ! 2 4 4 4 4 4 4
donde n = x1 + x2 + x3 + x4 y 0 ≤ θ ≤ 1.
[email protected] IEUV Estadística computacional 15-10-2014 6/25
 I II III IV El algoritmo EM

Estimación analítica

La función de log-verosimilitud es

l(θ) = log(L(θ|x)) = log(g(x|θ)),

l(θ) ∝ x1 log(2 + θ) + (x2 + x3 ) log(1 − θ) + x4 log(θ).
Al derivar l(θ) con respecto a θ se tiene

∂ x1 x2 + x3 x4
l(θ) = − + .
∂θ 2+θ 1−θ θ
∂
Resolviendo el sitema
∂θ l(θ) = 0, lo cual equivale a la ecuación
cuadrátical∗ (θ) = nθ2 + (2x2 + 2x3 − x1 + x4 )θ − 2x4 = 0 y
despejando θ se obtiene θ.
el EVM de
p
−(2x2 + 2x3 − x1 + x4 ) + (2x2 + 2x3 − x1 + x4 )2 + 8nx4
θbM V = .
2n
Por lo tanto, para x = (x√1 , x2 , x3 , x4 )T = (125, 18, 20, 34)T , una
15+ 53809
estimación es θbM V = 394 = 0, 626821498.
[email protected] IEUV Estadística computacional 15-10-2014 7/25
 I II III IV El algoritmo EM

Estimación utilizando la librería maxLik

l(θ) ∝ x1 log(2 + θ) + (x2 + x3 ) log(1 − θ) + x4 log(θ).

library(maxLik)
x<-c(125,18,20,34)
LogLikTheta<- function(param)
{
theta<-param[1]
x[1]*log(2+theta)+(x[2]+x[3])*log(1-theta)+
x[4]*log(theta)
}
EVM <- maxLik(logLik = LogLikTheta, start = c(theta=0.5))
#summary(mle)
coef(EVM)

theta
0.6268215
[email protected] IEUV Estadística computacional 15-10-2014 8/25
 I II III IV El algoritmo EM

Estimación mediante método de

Newton-Raphson

Recuerde que, si desea encontrar los valores de θ que satisfacen la

ecuación f (θ) = aθ2 + bθ + c = las puede obtener con el método
de N-R mediante el algoritmo iterativo

f (θ(j) ) a[θ(j) ]2 + bθ(j) + c

θ(j+1) = θ(j) − = θ (j) − .
f 0 (θ(j) ) 2aθ(j) + b
En nuestro ejemplo tendriamos

n[θ(j) ]2 + (2x2 + 2x3 − x1 + x4 )θ(j) − 2x4

θ(j+1) = θ(j) − .
2nθ(j) + (2x2 + 2x3 − x1 + x4 )
Reemplazando los datos

197[θ(j) ]2 − 15θ(j) − 68
θ(j+1) = θ(j) − .
394θ(j) − 15

[email protected] IEUV Estadística computacional 15-10-2014 9/25

 I II III IV El algoritmo EM

Estimación mediante el método de

Newton-Raphson

La EVM mediante al algoritmo de N-R es de θbN R = 0, 626821498.

Iteración j θ(j) δj = |θ(j) − θ(j−1) | l∗ (θj )

0 0,500000000 - -26,2500000000000
1 0,644230769 0,144230769 4,0980954142012
2 0,627071500 0,017159269 0,0580047808911
3 0,626821551 0,000249949 0,0000123075082
4 0,626821498 0,000000053 0,0000000000006
5 0,626821498 0,000000000 0,0000000000000

Por ejemplo,
197[θ(0) ]2 −15θ(0) −68 2 −15∗0,50−68
θ(1) = θ(0) − 394θ(0) −15 = 0, 50 − 197[0,50]
394∗0,50−15 = 0, 644230769
197[0,644230769]2 −15∗0,644230769−68
θ(2) = 0, 644230769 − 394∗0,644230769−15 = 0, 627071500

[email protected] IEUV Estadística computacional 15-10-2014 10/25

El algortitmo EM
 I II III IV El algoritmo EM

Primeros detalles del algoritmo EM

Se supone que X = (X1 , X2 , . . . , Xn ) VA iid con distribución

conjunta desde g(x|θ) y se quiere calcular

θb = arg máx L(θ|x).

Donde L(θ|x) = g(x|θ). Consideremos los datos completos w

provenientes una muestra aleatoria constituida por W = (X, Z),
donde W representa los datos completos, X los datos observados
y Z datos perdidos. La distribución conjunta de W es

f (w|θ) = f (x, z|θ) = k(z|θ, x)g(x|θ).

¾Cómo calculamos Lc (θ|w) = Lc (θ|x, z) si no conocemos z? Res-

puesta: No
c
conocemos z de L (θ|x, z), así que la supondremos
como variable aleatoria y calculamos una media.

[email protected] IEUV Estadística computacional 15-10-2014 12/25

 I II III IV El algoritmo EM

Primeros detalles del algoritmo EM

R
Considerando g(x|θ) = Z f (x, z|θ)dz, donde (X, Z) ∼ f (x, z|θ).
Entonces la distribución condicional de los datos perdidos z, dado
los datos observados x es

f (x, z|θ)
k(z|θ, x) = .
g(x|θ)

Además existe una relación entre la verosimilitud para los datos

completos Lc (θ|x, Z) y la verosimilitud para los datos observados
L(θ|x) dada por

Lc (θ|x, Z) = k(Z|θ, x)L(θ|x)

y la log-verosimilitud es

log Lc (θ|x, Z) = log k(Z|θ, x) + log g(x|θ)

[email protected] IEUV Estadística computacional 15-10-2014 13/25

 I II III IV El algoritmo EM

Primeros detalles del algoritmo EM

es decir,

log g(x|θ) = log L(θ|x) = log Lc (θ|x, Z) − log k(Z|θ, x)

Para un valor θ0 , calculando esperanza con respecto a k(Z|θ, x)

y utilizando la desigualdad de Jensen, se tiene

log L(θ|x) = Eθ0 [log Lc (θ|x, Z)] − Eθ0 [log k(Z|θ, x)] .
| {z } | {z } | {z }
Datos obs. Datos completos Datos perdidos

Al maximizar log L(θ|x) se debe ignorar el termino asociado solo

a los datos perdidos.

[email protected] IEUV Estadística computacional 15-10-2014 14/25

 I II III IV El algoritmo EM

Iteraciones

El valor esperado de la log-verosimilitud se denota por

Q(θ|θ0 , x) = Eθ0 [log Lc (θ|x, Z)].

El algoritmo EM comienza maximizando en cada iteración

Q(θ|θ0 , x).
Si θ
b(1) = arg máx Q(θ|θ0 , x), entonces θb(0) → θb(1) .
Se obtienen secuencias de estimadores {θb(j) }, donde

θb(j) = arg máx Q(θ|θb(j−1) , x).

Este esquema iterativo, en cada paso contiene un calculo de

esperanza y maximización.

[email protected] IEUV Estadística computacional 15-10-2014 15/25

 I II III IV El algoritmo EM

El algoritmo

Se comienza con una valor inicial θb(0) jado por el investigador.

Repita.
1 Calcule (paso E)

Q(θ|θb(m) , x) = Eθb [log Lc (θ|x, Z)],

(m)

donde la esperanza es con respecto a k(z|θb(m) , x) y establecer

m = 0.
2 Maximizar Q(θ|θb(m) , x) en θ y tomar (paso M)

θb(m+1) = arg máx Q(θ|θb(m) , x)

θ

y establecer m = m + 1.
Los parámetros que se encuentran en el paso M se usan para
comenzar el paso E siguiente, y así el proceso se repite. Es decir,
se ja el punto θb(m+1) = θb(m) .
[email protected] IEUV Estadística computacional 15-10-2014 16/25
Ejemplos del algortitmo EM
 I II III IV El algoritmo EM

Ejemplo: distribución de los datos observados

Recuerde que se presentaron datos en que n = 197 animales se dis-

tribuyen en cuatro categorías, de modo que el vector de datos ob-
servados de frecuencias es x = (x1 , x2 , x3 , x4 )T = (125, 18, 20, 34)T .
Se postula que provienen de una distribucion multinomial.

La función de distribución conjunta, para los datos observa-

dos, es:

  x1  x2 +x3  x4

n! 1 1 1 θ
g(x|θ) = + θ (1 − θ) .
x1 !x2 !x3 !x4 ! 2 4 4 4

El núcleo de la distribución de los datos observados es

g(x|θ) ∝ (2 + θ)x1 (1 − θ)x2 +x3 (θ)x4 .

[email protected] IEUV Estadística computacional 15-10-2014 18/25

 I II III IV El algoritmo EM

Ejemplo: distribución de los datos completos

Suponga que el dato con la mayor categoría proviene de dos ca-

tegorías z1 y z2 con probabilidades 1/2 y θ/4 , respectivamente.
De esta forma se introduce una variable latente (x1 = z1 + z2 ) re-
sultando cinco categorías dadas por w = (w1 , w2 , w3 , w4 , w5 )T =
(z1 , z2 , x2 , x3 , x4 )T , que representarían los datos completo. Ade-
más z = (z1 , z2 ) representa los datos perdidos.

La distribución conjunta de W es f (w|θ) = f (x, z|θ), donde

 z1  z2  x2 +x3  x4

n! 1 θ 1 θ
f (x, z|θ) = (1 − θ)
z1 !z2 !x2 !x3 !x4 ! 2 4 4 4

El núcleo de la distribución de los datos observados es

f (x, z|θ) ∝ (θ)z2 +x4 (1 − θ)x2 +x3 .

[email protected] IEUV Estadística computacional 15-10-2014 19/25

 I II III IV El algoritmo EM

Ejemplo: distribución condicional de los

datos perdidos

La distribución condicional de los datos perdidos z, dado los datos

observados x es
  z2  x1 −z2
f (x, z|θ) x1 θ θ
k(z|θ, x) = = 1− .
g(x|θ) y2 θ+2 θ+2
El nucleo de la distribución condicional de los datos perdidos z,
dado los datos observados x es

(θ)z2 +x4 (1 − θ)x2 +x3 θ z2

k(z|θ, x) ∝ = .
(2 + θ)x1 (1 − θ)x2 +x3 (θ)x4 (2 + θ)x1
En conclusión, la distribución de Z2 es binomial con n = x1 y
θ
p= θ+2 y por lo tanto

θx1
E(Z2 ) = .
θ+2
[email protected] IEUV Estadística computacional 15-10-2014 20/25
 I II III IV El algoritmo EM

Ejemplo: La etapa E

Primero la log-verosimilitud completa es

log Lc (θ|x, Z) = log f (x, Z|θ) = (Z2 +x4 ) log(θ)+(x2 +x3 ) log(1−θ).

La esperanza sería

Eb
θ(m)
[log Lc (θ|x, Z)] = Eθb [(Z2 +x4 ) log(θ)+(x2 +x3 ) log(1−θ)],
(m)

= Eθb [Z2 ] log(θ) + x4 log(θ) + (x2 + x3 ) log(1 − θ),

(m)

Por lo tanto,

θ(m) x1
Q(θ|θb(m) , x) = log(θ) + x4 log(θ) + (x2 + x3 ) log(1 − θ),
θ(m) + 2

[email protected] IEUV Estadística computacional 15-10-2014 21/25

 I II III IV El algoritmo EM

Ejemplo: la etapa M

Se deriva Q con respecto al parámetro θ y se iguala a cero

θ(m) x1
∂Q(θ|θb(m) , x) θ(m) +2 + x4 x2 + x3
= − = 0.
∂θ θ 1−θ
Resolviendo la ecuación (despejando θ)
θ(m) x1
θ(m) +2 + x4
θ= θ(m) x1
θ(m) +2 + x2 + x3 + x4

Remplazando los datos, x1 = 125, x2 = 18, x3 = 20 y x4 = 34,

resulta una formula recursiva para el EVM de θ mediante.
125θ(m)
θ(m) +2 + 34 159θ(m) + 68
θb(m+1) = 125θ(m)
= .
+ 72 197θ(m) + 144
θ(m) +2

[email protected] IEUV Estadística computacional 15-10-2014 22/25

 I II III IV El algoritmo EM

Ejemplo: la estimación

La EVM mediante el algoritmo EM es de θbEM = 0, 626821498.

Iteración j θ(j) δj = |θ(j) − θ(j−1) |

1 0,500000000 -
2 0,608247423 0,108247423
3 0,624321050 0,016073628
4 0,626488879 0,002167829
5 0,626777322 0,000288443
6 0,626815632 0,000038310
7 0,626820719 0,000005087
8 0,626821394 0,000000675
9 0,626821484 0,000000090
10 0,626821496 0,000000012
11 0,626821498 0,000000002
12 0,626821498 0,000000000

[email protected] IEUV Estadística computacional 15-10-2014 23/25

 I II III IV El algoritmo EM

Código R: Algoritmo EM Monte Carlo

emrao <- function(m,teta0,error) Iteración j θ(j)

{ 1 0.6080664

error1 <- 1 2 0.6243003

contador <- 1 3 0.6263641

while(error1 > error) 4 0.6269785

{ 5 0.6267242

tetaux <- teta0/(2+teta0) ... ...

149 0.6268228
zi <- rbinom(m,125,tetaux)
148 0.6271139
teta1 <- (mean(zi)+34)/(mean(zi)+72)
150 0.6267015
error1 <- abs((teta1-teta0)/teta1)
148 0.6268587
teta0 <- teta1
151 0.6268609
print(c(teta1,contador))
contador <- contador+1
La EVM mediante el
}
} algoritmo EMMC es

emrao(10000,0.5,0.00001) de θbM C = 0,6268609.

[email protected] IEUV Estadística computacional 15-10-2014 24/25

 Fin sesión 4
Gracias por asistir a esta presentación