Métodos Monte Carlo basados en cadenas de Markov

TFG: Métodos Monte Carlo basados en cadenas de Markov José Jiménez JIMENEZLUNA.COM Licensed under the Creative Commons Attribution-NonCommercial 4.0 Unported License (the “License”). You may not use this file except in compliance with the License. You may obtain a copy of the License at http://creativecommons.org/licenses/by-nc/4.0 . Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an “AS IS”BASIS,WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. Primera edición, Enero 2015 Abstract: Markov Chain Monte Carlo (or shortly MCMC) is a powerful method for sampling from high dimensional probability distributions. Here we present a swift introduction to the theory behind these methods as well as several applications of the Bayesian MCMC framework. These include, among others Bayesian Mixture Models, Bayesian Image Analysis or Text Mining. Implementations of solutions regarding these problems can be found in several programming languages. Índice general 1Motivación del problema .................................... 7 1.1 Problemática en inferencia bayesiana 7 1.2 Cálculo de esperanzas 8 2Conceptos previos .......................................... 9 2.1 Muestreo por rechazo 9 2.2 Cociente de uniformes 10 2.3 Integración Monte Carlo 10 2.4 Muestreo por importancia 11 2.5 Cadenas de Markov 12 2.6 Propiedades de las Cadenas de Markov 13 3Markov Chain Monte Carlo ................................. 15 3.1 Algoritmo de Metropolis-Hastings 15 3.1.1 Distribucionesproposición....................................... 16 3.2 Algoritmos de muestreo 17 3.2.1 AlgoritmodeMetropolis ........................................ 17 3.2.2 PaseoaleatorioMetropolis ...................................... 17 3.2.3 Muestreoconindependencia ................................... 17 3.2.4 Actualizandoenbloques ....................................... 18 3.3 Muestreo de Gibbs 19 3.4 Muestreo por rechazo adaptativo (ARS) 20 3.5 Otras consideraciones 21 3.5.1 Valoresiniciales ............................................... 21 3.5.2 Calentamiento ............................................... 21 3.5.3 Análisisdelasalida ............................................ 22 3.5.4 Tasadeconvergencia.......................................... 22 3.5.5 Estimacióndelavarianza ....................................... 22 4Modelos graficos y DAGs ................................... 25 4.1 Digrafos acíclicos 25 4.2 Grafo de independencia condicional 26 4.3 Ejemplos de aplicación 27 5Mixturas gaussianas ........................................ 33 5.1 Modelos de mixturas finitos 33 5.2 Usando MCMC 35 5.3 Dificultad en reasignación 39 5.4 Determinando el número de subpoblaciones 42 6Análisis de imagen ......................................... 45 6.1 Reconstrucción de imágenes 45 6.1.1 CamposaleatoriosdeMarkov ................................... 46 6.1.2 ModelodeIsing............................................... 46 6.1.3 ModelodePotts .............................................. 51 6.1.4 Sobre la constante de integración . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 6.1.5 Segmentacióndeimágenes..................................... 56 7Imputación múltiple ........................................ 63 7.1 Imputación usando ecuaciones encadenadas 63 7.1.1 Modelos de imputación univariante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 7.2 Probando el algoritmo 66 8Minería de texto ............................................ 69 8.1 Descubrimiento de temáticas 69 8.2 Asignación latente de Dirichlet 69 8.2.1 Formalización de la generación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 8.3 Inferencia y estimación usando muestreo de Gibbs colapsado 70 8.3.1 Formalización del aprendizaje . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 8.4 Ejemplos de aplicación 72 8.4.1 AnálisisdeperfilenTwitter ....................................... 73 8.4.2 Temática según autores clásicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 Bibliografía ................................................. 77 Problemática en inferencia bayesiana Cálculo de esperanzas 1. Motivación del problema Los métodos MCMC (Markov Chain Monte Carlo) surgen de la necesidad de simular el comportamiento de variables aleatorias y de estimar parámetros de las funciones de densidad/probabilidad de las mismas. El gran impulso a estas técnicas se las da mayormente (pero no únicamente) el enfoque estadístico bayesiano, donde la inferencia se realiza sobre lo que se denomina una función a posteriori, que denominaremos π(θ|x). Las siglas MCMC vienen marcadas por las cadenas de Markov y por la integración Monte Carlo. La inferencia bayesiana, en multitud de ocasiones necesita integrar sobre distribuciones de dimensión muy elevada (en muchas ocasiones, con cientos de parámetros). Existen métodos numéricos aproximados que producen buenas soluciones, pero que no escalan bien con la dimensión, siendo en muchos casos, computacionalmente intratables. De manera poco formal, la aplicación de técnicas MCMC consta de dos pasos: 1. Generar una muestra X1,...,Xn mediante una cadena de Markov cuya distribución estacionaria sea la buscada. 2. Tomar medias muestrales (integración Monte Carlo) y realizar inferencias sobre la muestra anteriormente citada Algunos términos no quedan explicados a estas alturas. No obstante, esta estructura general que quedará aclarada en los siguientes capítulos permite resolver muchos problemas. 1.1 Problemática en inferencia bayesiana Como hemos mencionado, en inferencia bayesiana, se realiza inferencia sobre una distribución a posteriori π(θ|x) . Bajo un enfoque bayesiano, no existe diferencia conceptual entre parámetros y valores observables, es decir, son en su totalidad cantidades aleatorias. Denominemos por x a los datos observados, y θ al conjunto de parámetros (nótese que pueden ser ambos multidimensionales). En estadística bayesiana, la función a posteriori puede expresarse de la siguiente manera, aplicando el teorema de Bayes: π(θ|x) = L(x|θ)π(θ) RL(x|θ)π(θ)(1.1) Donde L(x|θ) no es más que la función de verosimilitud de los datos asumiendo que siguen 8Capítulo 1. Motivación del problema una distribución parametrizada y π(θ) es lo que se denomina una distribución a priori. Ésta última expresa en inferencia bayesiana toda información previa o creencia que se tiene sobre el comportamiento de los parámetros del modelo. O lo que es lo mismo, en inferencia bayesiana los parámetros son también variables aleatorias, cuya información se incluye en un análisis a posteriori, tomando también como referencia la muestra xensayada. En muchísimas ocasiones (por no decir casi todas), esta distribución a posteriori se conoce únicamente salvo constante multiplicativa (que fuerza a que la distribución integre 1 sobre su dominio Ω). Por esta razón, generalmente se suele notar de la siguiente manera: π(θ|x)∝L(x|θ)π(θ)(1.2) Es decir, se dice que la función a posteriori es proporcional al numerador de 1.1. Una vez que tenemos determinada esta distribución, el enfoque bayesiano realiza inferencia sobre esperanzas de funciones de la misma. 1.2 Cálculo de esperanzas Como acabamos de mencionar, la problemática ahora surge de estimar E[f(x)] sobre la distribución a posteriori. En la mayoría de ocasiones trabajaremos sobre espacios paramétricos elevados, donde las soluciones analíticas pasan a ser un dolor de cabeza, y donde las numéricas no escalan bien. (De hecho, se pueden realizar simulaciones donde se comprueba que en espacios paramétricos con más de 20 dimensiones estos métodos suelen fallar estrepitósamente.) En el resto del trabajo seguimos una hoja de ruta encaminada a aprender cómo resolver estos problemas mencionados de acuerdo con los pasos citados en el punto anterior. Antes de ello, realizaremos un muy breve repaso a algunos métodos de simulación de variables aleatorias básicos (que serán útiles a lo largo del trabajo), luego recordaremos algunas propiedades básicas de las cadenas de Markov que nos darán soporte teórico sobre lo que queremos hacer, y para finalizar los conceptos previos hablaremos muy brevemente de integración Monte Carlo. Una vez finalizada la introducción, explicamos los métodos MCMC en general, dando varios algoritmos generales para resolución de dichos problemas (con ejemplos sencillos), algunos fundamentos teóricos para sustentar validación de resultados y algunos métodos de simulación más avanzados. La última parte (y la más pesada) del trabajo corresponde a varias aplicaciones prácticas de estos métodos. Algunos ejemplos son tratamiento de imágenes o mixturas. En esta parte del trabajo especialmente se hará un uso intensivo de lenguajes de programación, como pueden ser R o Python 2.7. Las librerías utilizadas en cada parte se detallan para posterior reproducibilidad. Todo el código del trabajo queda disponible en el repositorio github.com/hawk31/MCMC_tfg bajo licencia MIT. Muestreo por rechazo Cociente de uniformes Integración Monte Carlo Muestreo por importancia Cadenas de Markov Propiedades de las Cadenas de Markov 2. Conceptos previos En este capítulo repasamos algunas ideas previas para la mejor comprensión de las aplicaciones posteriores. Estas ideas básicas servirán para esbozar el marco de actuación bayesiano ante un problema determinado. 2.1 Muestreo por rechazo El muestreo por rechazo proporciona una manera muy eficiente de simular valores de una variable aleatoria f(x) . Suponemos en este sentido que no tenemos un algoritmo sencillo (como puede ser mediante inversión) para generar valores aleatorios de dicha distribución. La idea es utilizar una distribución instrumental g(x) , que acote absolutamente a f(x) , es decir f(x)<Mg(x) , con M>1 y de la que sepamos generar valores aleatorios de una manera fácil. El algoritmo funciona de la siguiente manera: Muestrear un valor xde g(x)yuperteneciente a una U(0,1). Si u<f(x) Mg(x) , aceptar x como valor aleatorio de f(x) . En caso contrario rechazar el valor y volver al primer paso. Nótese que este método supone que tanto f como g son evaluables y que se conoce la constante multiplicativa M hasta cierto punto. Este método, aunque eficiente, tiene algunas desventajas, como que la distribución g debe ser parecida en curtosis y simetría a f , de manera que en el segundo paso del algoritmo se acepte una proporción de veces aceptable. En caso de que esto no ocurra, desechamos muchas muestras. En este sentido, existen métodos un poco más complejos (Muestreo por rechazo adaptativo) que tratan este problema de una manera más atractiva. El muestreo por rechazo es precursor directo del algoritmo de Metrópolis-Hashtings, que además tiene la ventaja de escalar mucho mejor con la dimensionalidad de la simulación. Podemos realizar una pequeña simulación de valores aleatorios de procedentes de una N(0,1) utilizando como distribución instrumental una Cauchy(0,2), con constante M=3 . (Nota: si se desea ver la adecuación de esta distribución instrumental bastaría con representarla) 16 Capítulo 3. Markov Chain Monte Carlo La distribución q(.|.) puede tener cualquier forma y la cadena proporcionada por el algoritmo convergerá a la distribución estacionaria π(x) . Si bien esta condición es suficiente, la distribución proposición debe tener la misma dimensión que la estacionaria, y ser capaz de generar valores que se acepten. En otro caso la cadena puede pasar periodos largos de tiempo en un mismo estado. Nótese que en el momento que Xt ya pertenezca a la distribución estacionaria, todos los valores siguientes Xt+1,Xt+2,... pertenecerán igualmente a la misma distribución. Ejercicio 3.1 Para ilustrar cómo funciona el algoritmo de Metropolis-Hastings, utilizaremos un ejemplo de juguete. Supongamos que queremos muestrear valores de una distribución N(0,1) ,utilizando como distribución propuesta N(Xt,0,5) . Una posible implementación podría ser la siguiente. n = 1000 alea = numeric(n) alea[1]=0 #mu for (i in 2:n) {y = rnorm(1,alea[i-1],0.5) u = runif(1) alpha = min(1,(dnorm(y)*dnorm(alea[i-1],y,0.5))/(dnorm(alea[i-1])* dnorm(y,alea[i-1],0.5))) if(u<alpha) alea[i]=y else alea[i]=alea[i-1] } En la práctica, se suele realizar lo que se denomina ’traza’ de la simulación, que no es más que un gráfico de línea para comprobar la convergencia y el mezclado de la cadena. plot(alea,type="l")  3.1.1 Distribuciones proposición Como hemos notado antes, cualquier distribución proposición q(.|.) es suficiente para que la cadena converja a su distribución estacionaria. No obstante, la elección adecuada de esta distribución y su relación con la estacionaria garantizará una convergencia más rápida hacia esta 3.2 Algoritmos de muestreo 17 última. Es más, aún garantizada convergencia, la cadena podría ’mezclar’ de manera lenta (es decir, que se mueva lentamente por el soporte de la distribución objetivo). Por ello, escoger una distrubución propuesta adecuada se convierte en un problema. Existen varias formas canónicas que detallamos a continuación. 3.2 Algoritmos de muestreo Tradicionalmente, la literatura sobre MCMC ha hablado de muestreadores y algoritmos. No obstante, aunque aquí seguimos esas convenciones, conviene no olvidar que unos muestreadores no excluyen a otros. Como veremos un poco más adelante, es común combinarlos para construir una cadena de Markov con mejor rendimiento. Es por ello que tiene más sentido denominarlos como ’actualizadores MCMC’ más que muestreadores. Estos conceptos quedarán más asentados cuando expliquemos el muestreador de Gibbs. 3.2.1 Algoritmo de Metropolis El original propuesto por Metropolis en 1950 supone que utilizamos distribuciones proposición simétricas (esto es q(X|Y) = q(Y|X)). Con esta condicion 3.1, pasa a ser: α(Xt,Y) = min1,π(y) π(x)(3.2) 3.2.2 Paseo aleatorio Metropolis Si q(x,y) = f(y−x) para una densidad particular f , el algoritmo se denomina de paseo aleatorio Metropolis. Se suelen usar densidades f simétricas, para aceptar con probabilidades idénticas a las del apartado anterior. 3.2.3 Muestreo con independencia Supongamos a continuación que q(x,y) = f(y) , y por tanto los posibles candidatos a la cadena se generan independientemente del estado Xt de la cadena. En este caso, la probabilidad de aceptación puede escribirse como: α(x,y) = min1,w(y) w(x)(3.3) donde w(x) = π(x)/f(x) es lo que se denomina el peso de importancia. Este método guarda grandes similitudes con las ideas mostradas en el muestreo por importancia. La diferencia esencial entre los dos métodos es que el muestreador por importancia genera la masa de probabilidad sobre puntos con pesos altos, escogiéndolos frecuentemente. En contraste, el muestreo de independencia construye masa de probabilidad en puntos con altos pesos, permaneciendo en ellos largos periodos de tiempo. Ejercicio 3.2 Supongamos que queremos muestrear de una distribución Gamma de parámetros a,b cualesquiera usando el muestreador de independencia. Utilizaremos una distribución N(a/b,a/b2)como propuesta. Una posible implementación podría ser: #! /usr/bin/env RScript gammaSampler<-function (n, a, b) {mu <- a/b sig <- sqrt(a/b^2) 18 Capítulo 3. Markov Chain Monte Carlo vec <- numeric(n) x <- mu vec[1] <- x for (i in 2:n) { can <- rnorm(1, mu, sig) aprob <- min(1, (dgamma(can, a, b)/dgamma(x, a, b))/(dnorm(can, mu, sig)/dnorm(x, mu, sig))) u <- runif(1) if (u < aprob) vec[i] <- can else vec[i]<- vec[i-1] } vec } Nótese que hemos utilizado la media de la distribución objetivo como primer valor de la cadena para intentar conseguir una convergencia rápida. En la práctica esto no suele ser posible.  En la práctica, el muestreador de independencia puede funcionar muy bien o muy mal. Para que funcione bien, q(.) debería ser una buena aproximación a la distribución objetivo, aunque generalmente basta con que tenga colas más pesadas. En efecto, si q(.) no cumple esta condición, puede pasar mucho tiempo atascado en las colas de la distribución objetivo, proporcionando mal rendimiento. 3.2.4 Actualizando en bloques En los ejemplos que hemos estado realizando hasta ahora, hemos utilizado distribuciones con un número bajo de parámetros, para facilitar la comprensión. En la práctica desafortunadamente nos encontramos con situaciones en las que podemos llegar a tener cientos de parámetros. Es por ello que se han desarrollado procedimientos para actualizar cada uno de los componentes ’en bloque’ o secuencialmente. El método original propuesto por Metropolis se denomina singlecomponent Metropolis-Hastings. Supongamos que tenemos que actualizar X,1,....,X.h . Es decir, tenemos que actualizar h elementos, secuencialmente. Denotemos por X.−i al mismo vector sin el elemento i -ésimo. Para cada i de la iteración t+1 , actualizamos cada componente utilizando Metropolis-Hastings. El candidato Y.i se genera mediante distribución propuesta qi(Y.i|Xt.i,Xt.−i) . Esto quiere decir que la distribución propuesta puede depender del resto de componentes del vector y de sí misma en la iteración anterior. (Teniendo en cuenta que si se han actualizado algunos antes en la misma iteración t+1 pueden actualizarse los siguientes teniendo esto en cuenta). Por tanto, para cada componente del vector tendremos una probabilidad de aceptación al estilo anterior, con modificaciones pertinentes: α(X.−i,X.i,Y.i) = min1,π(Y.i|X.−i)qi(X.−i|Y.i,X.−i) π(X.i|X.−i)qi(Y.−i|X.i,X.−i)(3.4) En esta expresión π(.|X.−i) se denomina distribución totalmente condicionada. Esto no es más que la expresión fijando todos los parámetros menos el de interés, tratándolos como constantes. El uso de distribuciones totalmente condicionadas juega un papel muy importante en los métodos MCMC y en especial inferencia bayesiana, donde se aplican en modelos condicionales de 3.3 Muestreo de Gibbs 19 dependencia (ya los veremos en aplicaciones), y que en la mayoría de situaciones ofrecen una simplificación de 3.4. 3.3 Muestreo de Gibbs El muestreador de Gibbs cobra especial importancia justo después de este último apartado. El mismo nos proporciona una manera general y sencilla de actualizar los componentes del vector de parámetros mencionado utilizando distribuciones totalmente condicionadas. Utilizar Metropolis-Hastings para cada componente tiene la desventaja evidente de que tenemos que buscar h densidades propuesta, una para cada uno de los componentes. El muestreador de Gibbs propone lo siguiente: bajo las mismas condiciones que antes, la distribución propuesta de cada uno de los componentes será: qi(Y.i|X.i,X.−i) = π(Y.i|X.−i)(3.5) Es decir, se propone utilizar la distribución totalmente condicionada de un parámetro al resto de parámetros (valga la redundancia) como distribución propuesta. Esto, como hemos dicho antes pasa por fijar el resto de parámetros y tratarlos como constantes (valores iniciales o los de la iteración anterior). Finalmente, actualizamos la densidad resultante (univariante, esta vez) con el método que más convenga, ya sea Metropolis-Hastings u otro. Este método tiene gran aplicabilidad en todos los ejemplos prácticos que veremos al final del trabajo, ya que proporciona una manera directa de actualizar los componentes. Ejercicio 3.3 Para ejemplificar cómo funciona el muestreador de Gibbs, supongamos que queremos generar valores aleatorios de la siguiente distribución con 3 variables. f(y1,y2,y3)∝exp(−(y1+y2+y3+θ12y1y2+θ13y1y3+θ23y2y3)(3.6) Para aplicar el muestreador de Gibbs es necesario determinar las distribuciones totalmente condicionadas de y1,y2,y3. De esta manera: f(y1|y2,y3)∝exp(−y1(1+θ12y2+θ13y3)) (3.7) Nótese que podemos ignorar aquellos términos en los que no aparece y1 porque son constantes multiplicativas. Con un poco de ojo, y teniendo en cuenta que aquí y2 e y3 son constantes nos damos cuenta que y1|y2,y3∼Exp(1+θ12y2+θ13y3) . Por simetría y de manera idéntica, y2|y1,y3∼Exp(1+θ12y1+θ23y3) y y3|y1,y2∼Exp(1+θ12y1+θ23y2) . El ejemplo no estaría completo sin una simulación pertinente, para la cual supondremos unos θi j cualesquiera. niter = 10^4 theta = c(2,3,1/2) sim = matrix(NA,nrow=niter,ncol=3) sim[1,] = c(0.4,1.2,0.3) ## Valores iniciales (cualesquiera) for(i in 2:niter){ sim[i,1] = rexp(1,1+theta[1]*sim[i-1,2]+theta[2]*sim[i-1,3]) sim[i,2] = rexp(1,1+theta[1]*sim[i,1]+theta[3]*sim[i-1,3]) ## Nótese que aquí ya hemos actualizado el primer ## componente y lo usamos para acelerar convergencia. sim[i,3] = rexp(1,1+theta[2]*sim[i,1]+theta[3]*sim[i,2]) ## Idem 20 Capítulo 3. Markov Chain Monte Carlo } En este codigo hacemos uso de rexp , la función ya implementada en R para generar valores según una exponencial determinada. No obstante, no tenemos por qué tener un algoritmo eficiente en alguna de estas distribuciones condicionadas, por lo que podriamos usar Metropolis-Hastings o alguna de las técnicas ya descritas.  3.4 Muestreo por rechazo adaptativo (ARS) Aprovecharemos que ya hemos explicado el muestreador de Gibbs para explicar una técnica bastante usada en conjunto. Hasta ahora, cuando pretendíamos usar el muestreador de Gibbs suponíamos que teníamos que simular distribuciones totalmente condicionadas al resto, y que para ello podriamos recurrir a Metropolis-Hastings, al muestreo por rechazo o a cualquier otra técnica que proporcione buenos resultados. Si bien las anteriores técnicas funcionan bien, es necesario determinar en todas y cada una de ellas una densidad propuesta de la que generar los valores que luego serán aceptados o rechazados con una determinada probabilidad. El método que a continuación se propone es una técnica del tipo ’caja negra’ en el sentido de que no necesita esta determinación realizando una serie de hipótesis sobre la densidad de la que se pretende muestrear. Empezamos definiendo algunos conceptos. Supongamos que queremos muestrear de una densidad g(x) , la cual podemos conocer hasta una constante multiplicativa (constante de integración), continua y diferenciable en todo su dominio. Supongamos además que h(x) = ln g(x) es cóncava en todo D. Sea Tk= (x1,...,xk) , k puntos de abcisa en los que han sido evaluados tanto h(x) como h0(x) . A continuación definimos un contorno de rechazo en Tk como exp uk(x) , donde uk(x) es una función lineal definida a trozos definida por las tangentes de h(x) en Tk . Para j=1,...,k−1, las tangentes de xjyxj+1intersecan en: zj=h(xj+1)−h(xj)−xj+1h0(xj+1)+xjh0(xj) h0(xj)−h0(xj+1)(3.8) Esta función es una aproximación via contorno lineal por encima de g(x) .Por tanto, para x∈[zj−1,zj], definimos: uk(x) = h(xj)+(x−xj)h0(xj)(3.9) donde z0 es el extremo inferior de D (o −∞ si no está acotado inferiormente) y zk es el extremo superior de D (o ∞ si no está acotado superiormente.). Con estos conceptos, definimos la densidad: sk(x) = exp uk(x) RDexp uk(x0)dx0(3.10) Finalmente, definimos la función restricción inferior en Tk como explk(x) , donde lk(x) es una función lineal definida a trozos por debajo de g(x) , formada entre abcisas adyacentes en Tk . Para x∈[xj,xj+1]: lk(x) = (xj+1−x)h(xj)+(x−xj)h(xj+1) xj+1−xj (3.11) 3.5 Otras consideraciones 21 La concavidad de h(x) asegura que lk(X)≤h(x)≤uk(x) en todo el dominio. Tras todas estas definiciones, podemos poner en marcha el algoritmo. 1. Inicializar la abcisa en Tk . Si D no está acotado por la izquierda, escogemos x1 tal que h0(x1)>0 . Si D está no acotado por la derecha, entonces escogemos xk tal que h0(xk)<0 . 2. Calculamos las funciones uk(x),sk(x)ylk(x)para estos kpuntos. 3. Usamos sk(x)para muestrear un valor x∗ywde una uniforme unidad. 4. Si w≤exp(lk(x∗)−uk(x∗)) aceptar x∗. 5. En otro caso, si w≤exp(h(x∗)−uk(x∗)), aceptar x∗. En cualquier otro caso rechazarlo. 6. Si h(x∗) ha tenido que ser evaluado en este último test, incluímos x∗ en Tk para formar Tk+1y volvemos al segundo paso hasta que tengamos npuntos muestreados. Este algoritmo es el que se implementa en la mayoría de software especializado dedicado al muestreador de Gibbs, como pueden ser BUGS (o sus múltiples variantes) o JAGS , por lo que tiene gran importancia. En capítulos posteriores haremos uso de estas herramientas para algunas tareas. En R existe una implementación en el paquete ars . Ejercicio 3.4 Supongamos que queremos muestrear valores de una Beta(2,3) mediante Adaptive Rejection Sampling. Utilizamos el paquete mencionado. library(ars) #Muestrear de una Beta(2,3) n=20 f2<-function(x,a,b){(a-1)*log(x)+(b-1)*log(1-x)} f2prima<-function(x,a,b){(a-1)/x-(b-1)/(1-x)} mysample2<-ars(20,f2,f2prima,x=c(0.3,0.6),m=2, lb=TRUE,xlb=0,ub=TRUE,xub=1,a=2,b=3) mysample2  3.5 Otras consideraciones En este apartado consideraremos algunos aspectos a tener en cuenta cuando se están realizando simulaciones MCMC en general, desde cómo determinar unos valores iniciales adecuados hasta reglas para determinar un número de iteraciones adecuado. 3.5.1 Valores iniciales Si bien como hemos dicho los valores iniciales X0 no influyen en la distribución estacionaria de la cadena, si éstos no se escogen adecuadamente algunas cadenas pueden mezclar lentamente, teniendo que descartar muchas iteraciones como calentamiento. Una manera de determinar un buen valor inicial es correr varias cadenas simultáneamente con distintos valores iniciales (si se tiene la suficiente potencia computacional en paralelo). Si se dispone de información previa sobre la distribución estacionaria, un buen punto de partida podría ser la mediana. 3.5.2 Calentamiento El número de iteraciones de calentamiento, es decir, todas aquellas previas que se descartan por no pertenecer a la distribución de la que se pretende muestrear depende tanto de los valores iniciales como de la tasa de convergencia de la cadena (hablamos de esto último en el siguiente apartado). En mayor medida depende de cuán cerca estén las distribución propuesta y la estacionaria. La manera en la que se suele determinar el número m de iteraciones que se 22 Capítulo 3. Markov Chain Monte Carlo debe descartar se suele realizar mediante análisis gráfico de la traza de la cadena. Existen en la literatura algunos estimadores de m analíticos, pero no suelen verse con demasiada soltura en las aplicaciones probablemente porque el análisis gráfico proporciona la suficiente información. Otra cuestión relacionada es cuándo parar la cadena para realizar estimaciones Monte Carlo. La medida más usual para ello es correr como se ha sugerido antes varias cadenas en paralelo, de tal manera que si convergen rápidamente todas hacia una región determinada podemos asumir que ya ha alcanzado su distribución estacionaria y correr esta vez una cadena suficientemente larga para realizar estimaciones. Generalmente, las iteraciones pertenecientes al calentamiento se descartan en las estimaciones de medias y varianzas pues pueden sesgarlas. 3.5.3 Análisis de la salida Como se ha venido explicando mediante pinceladas hasta ahora, una salida de una simulación Monte Carlo proporciona suficiente información para: Realizar un gráfico de líneas de la traza para comprobar convergencia y calentamiento. Descartar estas últimas iteraciones en los siguientes análisis. Estimar medias y varianzas de la cadena. Podríamos incluso estimar intervalos de confianza, usando las estimaciones de varianza ahora mismo tomadas o de manera más sencilla, tomando los percentiles α 2 y 1−α 2 de la simulación. Rizando un poco el rizo, podriamos estimar las distribuciones marginales mediante alguna función núcleo K(X.i) . Una elección bastante corriente para esta función es la distribución totalmente condicionada de X.i|X.−i, como se explicó en el muestreador de Gibbs. 3.5.4 Tasa de convergencia Se dice que una cadena de Markov X es geométricamente ergódica si es aperiódica positiva recurrente y existe un λ∈(0,1)y una función V(.)tal que se cumple: ∑ j |Pi j −π(j)| ≤ V(i)λ(3.12) El λ más pequeño para el que exista una función V que satisfaga la condición anterior se denomina tasa de convergencia. Lo notaremos por λ∗ . Para entender mejor las implicaciones de las cadenas geométricamente ergódicas recurrimos al análisis espectral. Esto se escapa ampliamente del objetivo del trabajo, pero diremos que para las cadenas geométricamente ergódigas el principal autovalor es λ0=1 y el resto (finitos) están acotados por el círculo unidad. La velocidad a la que converge la cadena o λ∗ es por tanto dependiente del segundo autovalor más grande. 3.5.5 Estimación de la varianza Una de las consecuencias más importantes de las cadenas ergódicas es que permite la existencia de resultados del tipo teorema central del límite para medias ergódicas de la forma: fn−E[f(X)] →N(0,σ2)(3.13) Dicha convergencia es en distribución. Por tanto es de vital importancia estimar correctamente σ2 . Si bien pueden utilizarse estimadores simples como los propuestos en integración Monte Carlo, se han propuesto algunas alternativas más sofisticadas. 3.5 Otras consideraciones 23 Batch means La idea detrás de este procedimiento es correr una cadena de Markov en N=mn iteraciones, con nsuficientemente grande. Notemos por: Yk=1 n kn ∑ i=(k−1)n+1 f(Xi)(3.14) con k=1,...,m. De esta manera, una buena estimación de σ2podría venir dada por: ˆ σ2≈n m−1 m ∑ k=1Yk−fN(3.15) Ejercicio 3.5 Una posible implementación general de batch means podría ser la siguiente: n=1000 m=5 cadena = rexp(n*m,rate = 1/3) batch.means <- function(chain,f,m){ fN = mean(f(chain)) n = length(chain)/m li = numeric(m) for(i in 1:m){ li[i]=mean(f(chain[((i-1)*n+1):(i*n)])) } sigma = n/(m-1)*sum((li-fN)^2) return(sigma) } batch.means(cadena,identity,5)  Estimadores de ventana Otra opción para estimar σ2 es usar la función de autocovarianzas muestral. Sin embargo, este método produce estimadores inconsistentes según aumenta el retardo i , por lo que se propone una versión truncada del estimador, que puede expresarse de la siguiente manera: σ2≈ˆ γ0+2 ∞ ∑ i wn(i)ˆ γi(3.16) donde los wn(i) son pesos verificando |wn(i)|<1 y ∑iwn(i) = 1 . Una posible implementación rápida podria ser la siguiente: Ejercicio 3.6 cadena = rexp(5000,rate = 1/3) window.est <- function(cadena,weights){ aut = as.numeric(acf(cadena,plot=F,type="covariance")$acf) sigma = aut[1]+2*sum(pesos*aut[2:(length(pesos)+1)]) return(sigma) 24 Capítulo 3. Markov Chain Monte Carlo } pesos = rep(1/12,12) window.est(cadena,pesos)  Digrafos acíclicos Grafo de independencia condicional Ejemplos de aplicación 4. Modelos graficos y DAGs 4.1 Digrafos acíclicos Este capítulo será el más corto dedicado a aplicaciones. Los modelos gráficos se usan cada vez más frecuencia en modelos bayesianos en los que se aplica MCMC. Las relaciones entre las variables del modelo pueden representadas haciendo que los nodos en un grafo representen esas variables, y los ejes entre los nodos anteriores representando la relación (o ausencia) en términos de independencia condicional. Estos grafos, generalmente se presentan con una estructura jerárquica, con aquellos nodos que ejercen una influencia más directa sobre los datos colocados en la parte inferior, y según va disminuyendo dicha influencia colocando superiormente. Estos grafos proporcionan una manera fácil de interpretar la estructura condicional del modelo, simplificando la implementación del algoritmo MCMC en cuestión, indicando qué variable ejerce influencia sobre qué otra en su distribución. Estos grafos se denotan DAGs (Directed Acyclic Graphs). Consisten en una colección de nodos y ejes dirigidos, donde la dirección de estos ejes determina el sentido de la dependencia entre variables. Los nodos, a su vez, pueden ser de dos tipos: círculos si se trata de una variable desconocida y que por tanto hay que estimar simulando, o cuadrados si se trata de una variable conocida. Los DAGs se caracterizan por ser acíclicos, es decir, no podemos volver al mismo nodo de partida sea cual sea. Un ejemplo de DAG podría ser el de la figura 4.1. Modelos de mixturas finitos Usando MCMC Dificultad en reasignación Determinando el número de subpoblaciones 5. Mixturas gaussianas En este capítulo veremos una de las aplicaciones más interesantes que tienen las técnicas MCMC. Dentro del campo del aprendizaje no supervisado, podemos considerar las mixturas como parte de los algoritmos de clustering o conglomerados, donde dado un conjunto de datos el objetivo es describir si el mismo está compuesto de dos o más subpoblaciones bien diferenciadas. Concretamente, el enfoque que utilizaremos en este capítulo es intentar determinar si una muestra procede de dos o más distribuciones normales. Las mixturas son en realidad un caso particular de un conjunto de modelos denominados de variable latente o ausente. 5.1 Modelos de mixturas finitos En la mayoría de textos, una densidad de mixtura queda representada como una suma ponderada de funciones de densidad independientes (o una combinación convexa): k ∑ j=1 pjfj(x),;pj≥0; k ∑ j=1 pj=1 (5.1) Donde k es el número de subpoblaciones. En las situaciones más simples, las distribuciones fj se conocen y lo que interesa es estimar las cantidades pj , o en las probabilidades de que dado un elemento de la muestra, pertenezca a una subpoblación u otra. Generalmente, las distribuciones anteriores pertenecen a una distribución paramétrica, por lo que hay que incluir θ: k ∑ j=1 pjf(x|θj)(5.2) Dependiendo de la situación los objetivos de las mixturas pueden ser varios: puede ser estimar la pertenencia a los grupos de las observaciones z (clustering), para proporcionar estimadores de los parámetros de las distribuciones asociadas o incluso estimar el número de poblaciones subyacentes en el modelo. En cualquier caso, consideremos una muestra aleatoria simple x= (x1,...,xn) procedente de un modelo acorde a 5.1. Consideremos la función de verosimilitud asociada: 34 Capítulo 5. Mixturas gaussianas l(θ,p|x) = n ∏ i=1 k ∑ j=1 pjf(xi|θj)(5.3) A continuación empezamos el camino para intentar enfocar este problema a un marco de simulación MCMC. Para cualquier distribución a priori π(θ,p) , la distribución a posteriori de (θ,p)está disponible hasta constante multiplicativa (como de costumbre): π(θ,p|x)∝"n ∏ i=1 k ∑ j=1 pjf(xi|θj)#π(θ,p)(5.4) Para entender mejor cómo vamos a montar el modelo para la simulación, consideremos el enfoque de variable latente anteriormente mencionado. Para cada xi , diremos que tiene asociada una variable latente zi que indica la pertenencia de esa observación a una determinada subpoblación i . Una vez entendido esto, comenzamos definiendo algunas distribuciones condicionadas de manera adecuada: zi|p∼Mk(p1,..., pk)(5.5) En primer lugar, la variable latente zi condicionada a p se distribuirá según una distribución multinomial de parámetros definidos por p. xi|zi,θ∼f(.|θzi)(5.6) Las observaciones condicionadas a la subpoblación que pertenecen y a los parámetros de la subpoblación que pertenecen se distribuyen según f . En nuestro caso ésta función de densidad será la de una normal, pero el modelo es extensible a cualquier otra distribución (con mayor o menor dificultad). Con esta estructura, podemos definir de nuevo la función de verosimilitud: l(θ,p|x,z) = n ∏ i=1 pzif(xi|θzi)(5.7) Donde z= (z1,...,zn) .Por tanto la función a posteriori definida anteriormente pasaría a tener la siguiente forma, sustituyendo: π(θ,p|x,z)∝"n ∏ i=1 pzif(xi|θzi)#π(θ,p)(5.8) Ahora vamos a realizar una expansión conveniente para la expresión del modelo. Denotemos por Z=1,...,kn , el conjunto de los kn valores posibles del vector z . Z se puede descomponer trivialmente usando unión de conjuntos Z=∪τ j=1Zj . Siguiendo con este razonamiento, dado un vector de número de asignaciones a cada grupo (n1,...,nk) , definimos una serie de conjuntos de partición: Zj=(z: n ∑ i=1 Izi=1,..., n ∑ i=1 Izi=k)(5.9) Es un poco enrevesado seguir lo siguiente: lo anterior representa todas las posibles asignaciones dado un número de asignaciones (n1,...,nk) . Nombramos a cada uno de estos conjuntos de partición j=j(n1,...,nk) , usando el órden lexicográfico. Utilizando esto, la distribución a posteriori puede escribirse de forma cerrada: π(θ,p|x) = ∑ z∈Z π(θ,p|x,z) = τ ∑ i=1∑z∈Ziw(z)π(θ,p|x,z)(5.10) 5.2 Usando MCMC 35 Hemos enrevesado tanto la función a posteriori con dos finalidades: la primera encontrar una expresión cerrada para la función a posteriori, y lo segundo es para comprobar que w(z) representa la probabilidad marginal posterior de una asignación a z condicionada a x . Es por ello que un estimador Bayes de la distribución anterior es inmediato: E[θ,p|x] = τ ∑ i=1∑ z∈Zi w(z)E[θ,p|x,z](5.11) Esta descomposición tan poco natural que hemos desarrollado tiene bastante utilidad en otros métodos que preceden a las simulaciones MCMC que vamos a desarrollar a partir de ahora, concretamente al algoritmo EM. 5.2 Usando MCMC Llegados a este punto, nos planteamos cómo aplicar lo que hemos visto a las mixturas. La distribución totalmente condicionada de z|xestá disponible salvo constante multiplicativa: π(z|x,θ,p)∝ n ∏ i=1 pzif(xi|θzi)(5.12) A continuación falta determinar una serie de distribuciones conjugadas para las posteriori. Supongamos que p y θ son independientes a priori, entonces, dado z los vectores p y x son independientes: π(p|z,x)∝π(p)f(z|p)f(x|z)∝π(p|z)(5.13) De manera similar, θ es también independiente a posteriori de p|z y x , con densidad π(θ|z,x) . Aquí ya podemos aplicar el muestreo de Gibbs, simulando por una parte z condicionado sobre (p,θ) y sobre los datos (y al revés). De esta manera, nuestro algoritmo tendría que hacer lo siguiente: 1. Inicio: Escoger p(0)yθ(0)cualesquiera. 2. (t++) Para cada elemento de la muestra, generar z(t) i tal que P(zi=j|θ,p)∝p(t−1) jf(xi|θ(t−1) j) . 3. Generar p(t)según π(p|z(t)) 4. Generar θ(t)según π(θ|z(t),x) Para simular p , lo suyo es escoger una distribución conjugada en z . La complexidad de la simulación de los θj dependerá de la f con la que estemos trabajando (en nuestro caso, como son normales y existen distribuciones conjugadas para ambos parámetros no hay problema). Para lo primero que hemos mencionado, los zi se distribuyen como hemos dicho anteriormente según una multinomial M(p1,..., pk) , que permite encontrar fácilmente una distribución conjugada en z. De esta manera, la distribución a priori de p∼D(γ1,...,γk)será una Dirichlet, con densidad: Γ(γ1+...+γk) Γ(γ1)...Γ(γk)pγ1 1...pγk k(5.14) De esta manera, la distribución p|zserá Dirichlet también: p|z∼D(n1+γ1,...,nk+γk)(5.15) Ejercicio 5.1 R no trae una función de base para simular variables aleatorias según una distribución Dirichlet. Implementar una es fácil a través de valores aleatorios procedentes de distribuciones gamma o beta. Implementemos uno de manera rápida. 36 Capítulo 5. Mixturas gaussianas #!/usr/local/env RScript rdirichlet <- function(n=1,params=c(1,1)){ k = length(params) array = matrix(NA, nrow = n, ncol = k) for(i in 1:n){ support = rgamma(k,shape=params) array[i,] = support/sum(support) } return(array) }  Para las mixturas normales, las µj|σj,z,xson independientes, con distribución: µj|σj,z,x∼N εj(z),σ2 j nj+lj!(5.16) donde: εj(z) = ljεj+njxj(z) lj+nj (5.17) De manera similar: σ2 j|x,z∼IG(0,5(vj+nj),0,5sj(z)) (5.18) donde: sj(z) = s2 j+njs2 j(z)+ ljnj lj+nj (εj−xj(z))2(5.19) Después de tanta teoría, va siendo hora de que implementemos nuestro muestreo de Gibbs orientado a mixturas gaussianas con k componentes. En el código siguiente, asumimos que lj=1vj=20: Ejercicio 5.2 #!/usr/bin/env RScript # Gibbs Sampler (Gaussian mixtures) gibbsMixture <- function(data,k,max_iter=1000){ n = length(data) mu = mean(data) sig = var(data) # Preallocating data z = rep(NA,n) group_size = rep(NA, k) group_x = rep(NA, k) group_sum = rep(NA, k) mu_mat = sig_mat = pj_mat = matrix(NA, nrow=max_iter, ncol=k) mu_mat[1,] = rep(mu,k) 5.2 Usando MCMC 37 pj_mat[1,] = rep(1,k)/k sig_mat[1,] = rep(sig,k) ## Possible chunk for likelihood ####### ### Main loop for(i in 2:max_iter){ # z estimation for(j in 1:n){ p = pj_mat[i-1,]*dnorm(data[j],mu_mat[i-1,],sqrt(sig_mat[i-1,])) z[j] = sample(1:k, size = 1, prob = p) } # rest of parameters for(j in 1:k){ group_size[j] = length(which(z == j)) group_x[j] = sum(as.numeric(z == j)*data) group_sum[j] = sum(as.numeric(z==j)*(data-group_x[j]/group_size[j])^2) } mu_mat[i,] = rnorm(k, (mean(data) + group_x)/(group_size + 1), sqrt(sig_mat[i-1,]/group_size + 1)) sig_mat[i,] = 1/rgamma(k, .5*(20+group_size), var(data) + .5* group_sum + .5* group_size/(group_size + 1)*(mean(data) - group_x/group_size)^2) pj_mat[i,] = rdirichlet(params = group_size + 0.5) } class = setClass("MCMC mixture", slots=c(mu="matrix", sig="matrix", p="matrix", group="numeric")) res = class(mu = mu_mat, sig = sig_mat, p = pj_mat, group = z) return(res) }  38 Capítulo 5. Mixturas gaussianas Ejercicio 5.3 Usamos un marco de datos clásico como faithful, incluído en R para poner a prueba nuestro algoritmo. En primer lugar representamos los datos para hacernos una idea visual de si un modelo de mixturas es adecuado. data(faithful) hist(faithful$waiting, breaks = 20, col="lightblue", freq=F, main="Mixtura Gaussiana") Aparentemente parece que podemos modelar dicho conjunto de datos según una mixtura gaussiana de k=2 componentes. Aplicamos nuestro algoritmo y tomamos medias: set.seed(93) mixture = gibbsMixture(faithful$waiting,k = 2,max_iter = 1000) (mu = apply(mixture@mu, 2, mean)) [1] 55.11955 80.08556 (sig = apply(mixture@sig, 2, mean)) [1] 38.28160 33.79188 (p = apply(mixture@p, 2, mean)) [1] 0.367852 0.632148 Finalmente, podemos superponer sobre el histograma anterior la densidad de nuestra mixtura: curve(p[1]*dnorm(x,mu[1],sqrt(sig[1])),add=T, lwd=2) curve(p[2]*dnorm(x,mu[2],sqrt(sig[2])),add=T, col="red", lwd=2) 5.3 Dificultad en reasignación 39  Los resultados anteriores pueden llevar a una falsa sobrevaloración del modelo estadístico entre manos. Recordemos que el muestreador de Gibbs está condicionado a los valores iniciales que escojamos para la cadena. El condicionamiento sobre z implica que las cadenas y por tanto tanto θ como p son incapaces de realizar cambios drásticos sobre sí mismos y sobre las asignaciones en la siguiente iteración. Otras técnicas anteriores (como el algoritmo EM) también son sensibles a no realizar cambios drásticos sobre asignaciones y por tanto son igual de frágiles con respecto a este sentido. Una posible solución es utilizar el muestreo de Gibss en conjunto con alguna otra técnica MCMC, como Metropolis-Hastings. Podríamos usar en este sentido una probabilidad de aceptación con una adecuada distribución propuesta del tipo: π(θ0,p0|x) π(θ,p|x) q(θ,p|θ0,p0) q(θ0,p0|θ,p)(5.20) 5.3 Dificultad en reasignación Como se ha visto en la sección anterior, al muestreo de Gibbs le cuesta realizar cambios drásticos sobre reasignaciones en z . Una característica curiosa en los modelos de mixturas es que carece de un orden, es decir, dada una mixtura compuesta por dos densidades, es incorrecto decir que una de las dos es el primer componente de la misma o el segundo. De manera más técnica decimos que los parámetros θjno son marginalmente identificables. Por si fuera poco, en una mixtura de k componentes se puede probar que el número de modas en la verosimilitud es del orden de O(k!) . Esto es fácilmente comprobable ya que si (θ,p) es un máximo local en la verosimilitud, entonces una permutación de estos parámetros también sigue siéndolo. Es más, si usamos una distribución a priori sobre (θ,p) la cual es invariable bajo permutación de los índices, todas las distribuciones a posteriori marginales son idénticas, lo cual implicaría que los estimadores Bayes son los mismos. Este problema sobre reasignaciones tiene distintas soluciones. Por un lado, una solución 40 Capítulo 5. Mixturas gaussianas sencilla podría llevarse a cabo imponiendo una restricción de orden en los parámetros (por ejemplo, ordenando las medias). Esto se reduce en su totalidad a truncar la distribución a priori: π(θ,p)Iµ1≤...≤µk(5.21) Con esta reducción del espacio paramétrico pueden ocurrir sucesos inesperados, ya que las distribuciones a posteriori no tienen por qué respetar su topología, de tal manera que cuando ponemos simulaciones a funcionar los estimadores no tienen por qué acabar en una de las k! mencionadas anteriormente sino en una zona de silla de baja probabilidad. No obstante, si bien esta restricción puede suponer un cambio de rendimiento en las simulaciones, éstas no tienen por qué realizarse durante la misma sino después. En este sentido, si la restricción es sobre el orden en las medias, una vez que la simulación ha acabado, los componentes se pueden reasignar de acuerdo con este orden. Esto se puede realizar de la siguiente manera: dada una muestra de tamaño M de una simulación MCMC, definimos el estimador máximo a posteriori como: i∗=argmaxi=1,...,Mπ{(θ,p)(i)|x}(5.22) En palabras, es el valor simulado que maximiza la función de densidad a posteriori. Esta, como ya sabemos, no es necesario que incluya la constante de integración. Sin embargo, es bastante problable que este valor se encuentre en la vecindad de una de las k! posibles modas. Usaremos este valor óptimo (MAP) como pivote, en el sentido de que ordenaremos el resto de las otras iteraciones con respecto a dicha moda. En lugar de escoger este reorden según una distancia euclídea en el espacio paramétrico, la definimos en el espacio probabilístico de las asignacioens. Denotemos por Gk el conjunto de las kposibles permutaciones y τ∈Gk. Minimizamos en τuna distancia de entropía: h(i,τ) = n ∑ t=1 k ∑ j=1 P(zt=j|θ(i∗),p(i∗))×log P(zj=k|θ(i∗),p(i∗)) P(zt=j|τ(θ(i),p(i)))!(5.23) En definitiva, podemos definir nuestro algoritmo de reordenamiento pivotal como: En todas las iteraciones, calcular τi=argmin h(i,τ) (θ(i),p(i)) = τi{(θ(i),p(i))} Con este reordenamiento, en la mayoría de iteraciones las asignaciones quedan reasignadas a la misma moda, y por tanto eliminamos el problema de la identificación. Después de este reordenamiento, los estimadores Monte-Carlo siguen siendo los habituales. Realicemos la implementación de este reordenamiento en R: Ejercicio 5.4 Suponemos que tenemos una salida de la función gibbsMixture, la cual usamos en esta función: pivotalReor <- function(gibbsMix){ mu = gibbsMix@mu sig = gibbsMix@sig p = gibbsMix@p logpost = gibbsMix@logpost n = gibbsMix@n k = gibbsMix@k data = gibbsMix@data 5.3 Dificultad en reasignación 41 max_iter = gibbsMix@max_iter map_indices = order(logpost, decreasing = T)[1] map = list(mu = mu[map_indices,], sig= sig[map_indices,], p = p[map_indices,]) lili = matrix(NA, n, k) allocation = matrix(NA, n, k) for(t in 1:n){ lili[t,] = map$p*dnorm(data[t], mean = map$mu, sd=sqrt(map$sig)) lili[t,] = lili[t,]/sum(lili[t,]) } ordered_mu = matrix(NA, ncol=k, nrow=max_iter) ordered_sig = matrix(NA, ncol=k, nrow=max_iter) ordered_p = matrix(NA, ncol=k, nrow=max_iter) require(combinat) perma = permn(k) for(t in 1:1000){ entropies = rep(0, factorial(k)) for(j in 1:n){ allocation[j,] = p[t,]*dnorm(data[j], mean=mu[t,], sd=sqrt(sig[t,])) allocation[j,] = allocation[j,]/sum(allocation[j,]) for(i in 1:factorial(k)){ entropies[i] = entropies[i] + sum(lili[j,]*log(allocation[k,perma[[i]]])) } } best_ordering = order(entropies, decreasing=T)[1] ordered_mu[t,] = mu[t,perma[[best_ordering]]] ordered_sig[t,] = sig[t,perma[[best_ordering]]] ordered_p[t,] = p[t,perma[[best_ordering]]] } res = list(mu=ordered_mu, sig=ordered_sig, p=ordered_p) return(res) }  48 Capítulo 6. Análisis de imagen #include <RcppArmadillo.h> #include <RcppArmadilloExtensions/sample.h> #include <math.h> /* exp */ // [[Rcpp::depends(RcppArmadillo)]] using namespace Rcpp; // [[Rcpp::export]] int nei4(NumericMatrix x, int a, int b, int col){ int n = x.nrow(), m = x.ncol(); int nei = 0; a = a-1; b = b-1; if(a != 0){ if(x(a-1,b) == col){nei++;} } if(b != 0){ if(x(a,b-1) == col){nei++;} } if(a != (n-1)){ if(x(a+1,b) == col){nei++;} } if(b != (m-1)){ if(x(a,b+1) == col){nei++;} } return nei;} // [[Rcpp::export]] NumericMatrix isingSampler(NumericMatrix x, int max_iter, double beta){ int n = x.nrow(), m = x.ncol(); int i; int k; int l; NumericVector rows(n); rows = Rcpp::seq(1,n); NumericVector cols(m); cols = Rcpp::seq(1,m); NumericVector pos(2); pos = Rcpp::seq(0,1); NumericVector permut_rows(n); NumericVector permut_cols(m); NumericVector prob(2, 0.5); int n0; int n1; 6.1 Reconstrucción de imágenes 49 int pos1; int pos2; double res; for(i = 0; i < max_iter; i++){ permut_rows = RcppArmadillo::sample(rows,n,0); permut_cols = RcppArmadillo::sample(cols,m,0); for(k = 1; k <= n; k++){ for(l = 1; l <= m; l++){ n0 = nei4(x, permut_rows[k-1], permut_cols[l-1], 0); n1 = nei4(x, permut_rows[k-1], permut_cols[l-1], 1); prob[0] = exp(beta*n0); prob[1] = exp(beta*n1); pos1 = permut_rows[k-1]-1; pos2 = permut_cols[l-1]-1; double suma = prob[0] + prob[1]; prob[0] = prob[0]/suma; prob[1] = prob[1]/suma; res = RcppArmadillo::sample(pos, 1, 0, prob)[0]; //res = Rcpp::Function(sample(pos,1,0,prob)); Rcpp::Function(gc()); x(pos1, pos2) = res; } } Rcpp::Function(gc()); } return x; } Un asunto de interés es comprobar cómo se comporta el algoritmo dependiendo del valor de βescogido. Podemos realizar una prueba empírica rápida: betas = seq(0.2, 1.8, by = 0.3) 50 Capítulo 6. Análisis de imagen par(mfrow=c(3,2)) set.seed(39) x = matrix(sample(c(0,1),50*50,rep=T), nrow=50) for(beta in betas){ set.seed(93) y = isingSampler(x, 1000, beta) image(y, main=paste("beta=",beta)) } Que produce el siguiente resultado: 10 20 30 40 50 10 20 30 40 50 beta= 0.2 1:50 1:50 10 20 30 40 50 10 20 30 40 50 beta= 0.5 1:50 1:50 10 20 30 40 50 10 20 30 40 50 beta= 0.8 1:50 1:50 10 20 30 40 50 10 20 30 40 50 beta= 1.1 1:50 1:50 10 20 30 40 50 10 20 30 40 50 beta= 1.4 1:50 1:50 10 20 30 40 50 10 20 30 40 50 beta= 1.7 1:50 1:50 Como se puede comprobar, cuanto mayor es β , el modelo tiende a concentrarse bajo distribuciones de un solo color, y por tanto más homogénea la imagen.  6.1 Reconstrucción de imágenes 51 6.1.3 Modelo de Potts El modelo de Potts es una generalización natural al modelo de Ising cuando la imagen tiene más de dos colores, por ejemplo G . Esta vez notamos por ni,g al número de vecinos de i con color g , es decir, ni,g=∑j∼iIxj=g . La distribución totalmente condicionada de xi se escoge de manera idéntica como: π(xi=g|x−i)∝exp(βni,g)(6.5) La densidad conjunta del modelo de Potts, como ya esperamos tiene la forma: π(x)∝exp β∑ j∼i Ixj=xi!(6.6) Aquí nos encontramos con el mismo problema que en la sección anterior, y es que no conocemos el valor de β , y por tanto no podemos calcular la función de verosimilitud. Por otra parte, con un β grande, tal y como pasaba antes, es difícil que un valor se actualice y se reduce la velocidad de convergencia. Se propone una modificación basada en pasos Metropolis-Hastings que fuerza actualizaciones en cada paso. 1. Inicialización: Para cada i∈I, generar x(0) i∼UD(1,...,G) 2. Iteración t>1. Generar u, una permutación aleatoria de los elementos de I. 3. Para ldesde 1 hasta |I|: generar xu,l∼UD1,2,...,x(t−1) u,l−1,x(t−1) u,l+1,...,G generar número de vecinos n(t) ul,gy generar probabilidad de aceptación pl=maxexp(βnul,xul )/exp(βnul,xul )(t),1 si se acepta, sustituir xul por x(t) ul Ejercicio 6.2 De nuevo seguimos con la misma mecánica aplicada en el modelo anterior. Una implementación ineficiente del algoritmo en R podría ser la siguiente: pottsSampler <- function(x, num_col=2, max_iter=1000, beta){ n = dim(x)[1]; m = dim(x)[2] for(i in 1:max_iter){ permut = sample(1:(n*m)) cat("Iteracion", i, "\n") for(k in 1:(n*m)){ xcur = x[permut[k]] a = (permut[k]-1)%%n + 1 b = (permut[k]-1)%/%n + 1 xtilde = sample((1:num_col)[-xcur],1) prob = beta*(nei4(x,a,b,xtilde)-nei4(x,a,b,xcur)) if(log(runif(1))<prob){ x[permut[k]] = xtilde } } 52 Capítulo 6. Análisis de imagen } return(x) } Y una varios órdenes más eficiente, en C++ . #include <RcppArmadillo.h> #include <RcppArmadilloExtensions/sample.h> #include <math.h> // [[Rcpp::depends(RcppArmadillo)]] using namespace Rcpp; // [[Rcpp::export]] NumericMatrix pottsSampler(NumericMatrix x, int num_col, int max_iter, double beta){ int n = x.nrow(), m = x.ncol(); int i; int xcur; int k, j; int a, b; int xtilde; double prob; NumericVector pixels(n*m); pixels = Rcpp::seq(1, n*m); for(i = 0; i < max_iter; i++){ NumericVector permut(n*m); permut = RcppArmadillo::sample(pixels,n*m,1); for(k = 0; k < (n*m); k++){ NumericVector colors(num_col); colors = Rcpp::seq(1, num_col); xcur = x[permut[k] - 1]; a = remainder((permut[k]-1),n) + 1; b = (permut[k]-1)/n + 1; colors.erase(xcur - 1); xtilde = RcppArmadillo::sample(colors, 1, 0)[0]; prob = beta*(nei4(x, a, b, xtilde) - nei4(x, a, b, xcur)); double alea = runif(1)[0]; if(log(alea) < prob){ x[permut[k] - 1] = xtilde; } } 6.1 Reconstrucción de imágenes 53 } return x; } De manera similar a lo que hicimos para el modelo de Ising, podemos comprobar de manera empírica el comportamiendo del modelo teniendo en cuenta el valor de β. betas = seq(0.2, 1.4, by = 0.4) par(mfrow=c(2,2)) set.seed(39) x = matrix(sample(1:4,50*50,rep=T), nrow=50) for(beta in betas){ set.seed(93) y = pottsSampler(x, num_col=4, 1000, beta) image(x=1:50, y=1:50,z=y, main=paste("beta=",beta)) } 10 20 30 40 50 10 30 50 beta= 0.2 1:50 1:50 10 20 30 40 50 10 30 50 beta= 0.6 1:50 1:50 10 20 30 40 50 10 30 50 beta= 1 1:50 1:50 10 20 30 40 50 10 30 50 beta= 1.4 1:50 1:50 De manera similar, cuanto mayor es el parámetro, más tendencia se tiene hacia imágenes más homogéneas.  6.1.4 Sobre la constante de integración En las secciones anteriores hemos supuesto que β , y que por tanto Z(β) eran conocidos. Como cabe esperar, esto no se suele producir en la mayoría de las ocasiones. El manejo de esa constante de integración ha dado lugar a mucha literatura, pues es un problema difícil. Aquí veremos una metodología relativamente sencilla para averiguarla, denominada muestreo por caminos. Esta técnica está basada en una representación de la derivada de dicha constante: 54 Capítulo 6. Análisis de imagen dZ(β) dβ=∑ x S(x)exp(βS(x)) (6.7) Esta derivada puede expresarse de manera conveniente como una esperanza: dZ(β) dβ=Z(β)∑ x S(x)exp(βS(x) Z(β)=Z(β)E[S(x)] (6.8) y por tanto, tomando logaritmo antes: dlogZ(β) dβ=E[S(x)] (6.9) Esta representación permite que podamos expresar el ratio Z(β1 Z(β0)como una integral: log(Z(β1)/(Z(β0)) = Zβ1 β0 E[S(x)]dβ(6.10) Esta última ecuación es lo que se denomina identidad del muestreo por caminos. Con esta expresión, podemos recurrir a procedimientos de simulación estándares para su aproximación. La integral puede aproximarse por métodos numéricos, y para un valor de β , E[S(x)] = f(β) puede simularse usando el modelo de Potts, por ejemplo. Finalmente aproximamos f(β) por una función lineal a trozos para integrar fácilmente. Ejercicio 6.3 Realicemos una aproximación de f(β) usando R. Para un valor determinado de β: pathSampling <- function(x, ncol=2, max_iter = 10^3, beta){ n = dim(x)[1]; m = dim(x)[2] S = 0 for(i in seq_len(max_iter)){ if(i%%100==0){ cat("Iteración ",i,"\n") } s = 0 rows = sample(1:n) cols = sample(1:m) for(k in seq_len(n)){ for(l in seq_len(m)){ n0 = nei4(x, rows[k], cols[l], x[rows[k], cols[l]]) col = sample(1:ncol, 1) n1 = nei4(x, rows[k], cols[l], col) if(log(runif(1))< (beta*(n1-n0))){ x[rows[k], cols[l]] = col n0 = n1 } s = s+n0 } } if(2*i > max_iter){ 6.1 Reconstrucción de imágenes 55 S=S+s } } return(2*S/max_iter) } Y para generar una secuencia, en intervalos tan pequeños como se quiera: generatefbeta <- function(x, max_iter=10^3, ncol=2){ search = seq(0.1, 2, by=0.1) Z = seq_along(search) i = 1 for(beta in search){ cat("Beta: ", beta, "\n") Z[i] = pathSampling(x, ncol, max_iter, beta) i = i+1 } plot(search, Z, main="f(beta) approx.", type="l") return(Z) } El problema en eficiencia es que para imágenes de juguete nuestro algoritmo está bien, pero como venimos acostumbrando en este capítulo, sería buena idea implementar pathSampling en C++ . double pathSampling(NumericMatrix x, int ncol, int max_iter, double beta){ int n = x.nrow(), m = x.ncol(); double S = 0.0, s = 0.0; int i, k, l, n0, n1, col; NumericVector rows(n); rows = Rcpp::seq(1, n); NumericVector cols(n); cols = Rcpp::seq(1, m); NumericVector colors(ncol); colors = Rcpp::seq(1, ncol); for(i = 0; i < max_iter; i++){ s = 0.0; rows = RcppArmadillo::sample(rows, n, 0); cols = RcppArmadillo::sample(cols, m, 0); for(k = 0; k < n; k++){ 56 Capítulo 6. Análisis de imagen for(l = 0; l < n; l++){ n0 = nei4(x, rows[k], cols[l], x(rows[k]-1, cols[l]-1)); col = RcppArmadillo::sample(colors, 1, 0)[0]; n1 = nei4(x, rows[k], cols[l], col); if(log(runif(1)[0]) < (beta*(n1-n0))){ x(rows[k]-1, cols[l]-1) = col; n0 = n1; } s = s + n0; } } if(2*i > max_iter){ S = S + s; } } return 2*S/max_iter; }  6.1.5 Segmentación de imágenes Después de todas estas secciones ya tenemos todos los materiales disponibles para la tarea que queríamos realizar. Consideraremos las imágenes como objetos estadísticos, pero no como tal, sino que consideraremos que tenemos una imagen con cierta distorsión y , es decir que el color de gris (o cualquier otra escala, por conveniencia) se presenta con cierta perturbación. El objetivo de la segmentación de imágenes es, dada esta imagen distorsionada, realizar un procedimiento de ’clustering’ de píxeles atendiendo únicamente a la estructura de dependencia espacial de la estructura. Esta dependencia espacial ’real’ de píxeles se nota por x , donde generalmente y puede tomar valores reales positivos y x solo discretos, por conveniencia. Estamos pues interesados en conocer la distribución a posteriori de x dada y , es decir π(x|y)∝f(y|x)π(x) . En esta expresión f(y|x) representa la verosimilitud de los datos y la función link entre la imagen original y su clasificación, mientras que π(x) representa nuestra información a priori (o propiedades deseadas) sobre el comportamiento de la imagen ’real’. Generalmente esta distribución a priori suele ser el modelo de Potts ya estudiado con G categorías: π(x|β) = 1 Z(β)exp β∑ i∼j Ixj=xi!(6.11) Dado x , y por facilidad computacional, asumimos que las observaciones en y son variables independientes normales. Esto se hace porque es más fácil parametrizar que una distribución multinomial que tome valores en 0,...,255: 6.1 Reconstrucción de imágenes 57 f(y|x,σ2,µ1,...,µG) = ∏ i∈I 1 (2πσ2)(1/2)exp−1 2σ2(yi−µxi)(6.12) Para los parámetros de estas normales se suelen utilizar distribuciones a priori uniformes: β∼U(0,2)(6.13) µ∼U(0≤µ1≤... ≤µG≤255)(6.14) π(σ2)∝σ−2I[0,∞)(σ2)(6.15) Esta última no es más que una distribución uniforme en el logaritmo de σ . Los µi no tienen por qué ordenarse, pero así evitan el problema del reordenamiento pivotal explicado en el capítulo de las mixturas. La distribución conjunta a posteriori es la siguiente: π(x,β,σ2,µ|y)∝π(β,σ2,µ)×1 Z(β)exp β∑ j∼i Ixj=xi! ×∏ i∈I 1 (2πσ2)1/2exp−1 2σ2(yi−µxi)2(6.16) Comenzamos a continuación a construír las diversas distribuciones totalmente condicionadas para el trabajo del muestreo de Gibbs. La de xi: P(xi=g|y,β,σ2,µ)∝exp β∑ i∼j Ixj=g−1 2σ2(yi−µg)2!(6.17) Aquí puede comprobarse que, una vez que x es conocido, los elementos de distinta categoría se separan y el resto de parámetros pueden simularse independientemente condicionados a x,y y σ2 . Si notamos por ng=∑i∈IIxi=g y sg=∑i∈IIxi=gyi , la distribución totalmente condicionada de µg es una distribución normal truncada en [µg−1,µg+1] (por razones obvias, µ0=0,µg+1=255 ), de media sg/ng y varianza σ2/ng . La distribución condicionada de σ2 es una gamma inversa con parámetros |I|2/2 y ∑i∈I(yi−µxi)2/2. Finalmente, la distribución totalmente condicionada de βes aquella que: π(β|y)∝1 Z(β)exp β∑ j∼i Ixj=xi!(6.18) . Ejercicio 6.4 Una implementación rápida en R de la segmentación de imágenes bayesiana podría ser la siguiente: bay <- function (y, ncol=2, max_iter=10^3) {numb = dim(y)[1] x=0*y mu = matrix(0, max_iter, 6) sigma2 = rep(0, max_iter) mu[1, ] = c(35, 50, 65, 84, 92, 120) sigma2[1] = 20 beta = rep(1, max_iter) 64 Capítulo 7. Imputación múltiple Estos son solo algunos de los problemas que podrían surgir en nuestro análisis. Anteriormente, existía también un método de imputación múltiple basado en MCMC que los datos venían definidos por algún modelo de densidad multivariante P(Y|θ) . La dificultad de este método es evidente en el momento que tenemos un marco de datos con variables de distinto tipo. Aunar todas esas variables bajo un mismo modelo teórico multivariante (como por ejemplo hace el algoritmo EM y la normal multivariante) es por un lado irreal y por otro difícil. FCS es un método más natural en el sentido de que no partimos de una densidad multivariante. En su lugar la definimos implícitamente especificando una densidad univariante para cada variable del estilo P(yj|y−j,θj) . Esta densidad es la usamos para imputar yaus j dado y−j mediante algún tipo de regresión (generalmente lineal o logística) sobre los casos yobs j . Se realizan tantos paseos sobre todas las densidades como iteraciones del algoritmo sean necesarias. Las ventajas de FCS sobre otros métodos son evidentes: evitar tener que especificar directamente una distribución multivariante permite mucha flexibilidad a la hora de trabajar los datos. Es decir, convertimos un problema k dimensional en k problemas unidimensionales. Por otra parte, la idea de especificar un modelo de imputación diferente para cada variable es bastante natural. No obstante, FCS también cuenta con algunas dificultades. El tener que especificar un modelo diferente y adecuado para cada variable puede ser por un lado pesado para el investigador, y por otro es computacionalmente mucho más exigente. Por otra parte, evaluar la calidad de las imputaciones puede resultar difícil ya que la densidad conjunta teórica no tiene por qué existir. En cualquier caso, comencemos a definir el método. Supongamos que tenemos Y= (Y1,...,Yk) un vector de k variables aleatorias con una distribución P(Y|θ) . Asumimos que la anterior distribución queda totalmente definida por θ . El procedimiento general comprendería los siguientes pasos: 1. Determinar la distribución a posteriori p(θ|yobs)de θusando los datos observados yobs. 2. Muestrear un valor θ∗de p(θ|yobs) 3. Por último muestrear un valor y∗ de p(yaus|yobs,θ=θ∗) , la distribución condicional a posteriori de yaus dado θ∗. Generalmente este procedimiento es pesado y solamente se repiten los pasos 2 y 3 un número modesto de veces (generalmente suelen ser suficientes entre 5 y 10 iteraciones). Lógicamente desde el marco teórico univariante esto parece fácil, pero si Y es multivariante (caso k>1 ) lo lógico es utilizar un marco de simulación basado en el muestreo de Gibbs. El paso realmente complicado es el primero, donde necesitamos determinar la distribución a posteriori p(θ|yobs) . En nuestro caso, aplicamos el muestreo de Gibbs muestreando de distribuciones condicionales de la forma: P(Y1|Y−1,θ1)(7.1) ... (7.2) P(Yk|Y−k,θk)(7.3) Es decir, utilizamos los parámetros θi solo en las densidades univariantes donde intervienen respectivamente. En otras palabras, estos parámetros no serían necesariamente el producto de una distribución conjunta P(Y|θ) . De manera más clara, una iteración del muestreo de Gibbs comprendería los siguientes pasos: 7.1 Imputación usando ecuaciones encadenadas 65 θ∗ 1∼P(θ1|yobs 1,yt−1 2,...,yt−1 k)(7.4) y∗(t) 1∼P(yaus 1|yobs 1,yt−1 2,...,yt−1 k,θ∗(t) 1)(7.5) ... (7.6) θ∗ k∼P(θk|yobs k,yt 2,...,yt k−1)(7.7) y∗(t) k∼P(yaus k|yobs k,yt 2,...,yt k−1,θ∗(t) k)(7.8) (7.9) Y si nos fijamos, en realidad ya podemos recurrir a modelos de imputación donde y es univariante para muestrear de las anteriores distribuciones. En la siguiente sección se detallan varios modelos dependiendo de la escala de la variable, con lo que el problema quedaría completamente solucionado. En ningún caso se utiliza información (cuál se iba a usar, en cualquier caso) de yaus j para muestrear de θ∗(t) j. Las simulaciones del modelo pueden ser paralelizadas fácilmente. 7.1.1 Modelos de imputación univariante Para los pasos del muestreo de Gibbs es necesario recurrir a modelos de imputación univariante. Dependiendo de la distribución de y|x, los algoritmos son los siguientes: Datos normalmente distribuídos (escala) Para este caso, recurrimos a un modelo de regresión lineal ordinario. Suponemos media βx , varianza σ2 como parámetros. Por otra parte descomponemos x= (xobs,xmis) como la matriz n×k de datos y n=nobs +naus para la variable a predecir. El algoritmo consta de los siguientes pasos: 1. Estimar β por una regresión ordinaria lineal de y sobre x . Concretamente por: b= (xobs0xobs)−1xobs0yobs 2. Muestrear g∼χ2(nobs −p) 3. Estimar σ2∗= (yobs −xobsb)0(yobs −xobsb/g) 4. Muestrear w1∼N(0,Ik) 5. Calcular b∗=b+σ2∗w1V1/2 , donde V1/2 es la matriz triangular superior de una decomposición Cholesky de V= (xobs0xobs)−1 6. Muestrear w2∼N(0,Iaus n) , donde este último término representa una matriz triangular de tamaño ndonde los elementos donde yiesté ausente valen 1 y 0 en caso contrario. 7. Imputar y∗=xausb∗+w2σ∗ Este algoritmo debería ser bastante robusto a falta de normalidad, pero no obstante se proponen un par de alternativas para tener en cuenta este fenómeno: Predictive mean matching . Basta sustituir el último paso calculando yaus =xausb∗ . Para cada valor ausente, imputar ese valor por el valor de yobs =xobsb∗ que más cercano se encuentre ayaus. Hot-deck Reemplazar el penúltimo paso el muestreo por uno de naus valores con reemplazamiento del conjunto de nobs residuos estandarizados. 66 Capítulo 7. Imputación múltiple Datos de respuesta binaria Recurrimos a un modelo de regresión logística. Asumimos P(yi|β,xi) = exiβ (1+exiβ)yi 1−exiβ (1+exiβ)1−yi! . A continuación, las imputaciones se obtienen de la siguiente manera: 1. Obtener por un método numérico habitual (Newton-Raphson por ejemplo) un estimador b de βy de Var[β](un hessiano) usando solo los datos completos. 2. Muestrear b∗∼N(b,Var[b]) 3. Para cada observación ausente, calcular wi=exib∗ 1+exib∗ 4. Para cada observación ausente, muestreamos un ui∼U(0,1) . Si ui>wi , imputamos esa observación como yi=1, en otro caso, yi=0. Datos categóricos Este caso es el más sencillo porque no es más que una generalización natural del anterior. Notemos las categorías por 0,...,s−1 . La distribución de y puede caracterizarse por log(P(y= j|x)/P(y=0|x)) = βjx . En este sentido, el modelo para y no es más que una serie de regresiones logísticas apiladas considerando una categoría base (one vs rest). 1. De manera similar, utilizar un modelo de regresión logística multinomial, estimar b=ˆ β y Var(b) = ˆ Var(β). 2. Muestrear b∗∼N(b,Var(b)) 3. Para cada observación ausente, calcular πaus i,j=e−b∗ jxi 1+∑s−1 v=1eb∗ vxi 4. Para cada observación ausente xi, muestrear de 0,...,s−1 con probabilidad πaus i,j. 7.2 Probando el algoritmo En esta corta sección nos dedicaremos a probar con varios marcos de datos de distinta naturaleza cómo se comporta el algoritmo a nivel de calidad de imputaciones. Ejercicio 7.1 En primer lugar probaremos con el marco de datos mtcars del paquete datasets . Este marco se caracteriza por tener variables de todo tipo, por lo que se presta muy bien a analizar los errores cometidos. Por otra parte, es un marco de datos pequeños, por lo que servirá para evaluar la calidad de las imputaciones cuando hay pocas observaciones involucradas. library(mice) data = mtcars data$am = factor(data$am) data$vs = factor(data$vs) datosoriginales = data Introducimos de manera aleatoria datos ausentes: n = nrow(data) m = ncol(data) p = 0.15 7.2 Probando el algoritmo 67 ## Esperamos aprox n*m*p datos ausentes for(i in 1:n){ for(j in 1:m){ u = runif(1) if(u < p){ data[i,j] = NA } } } expected = n*m*p length(which(is.na(data))) Corremos el algoritmo, usando para todas las variables Predictive Mean Matching: imputation = mice(data, m=10, maxit=100, method = rep("fastpmm", 11)) completos = complete(imputation, action = 5) Finalmente, evaluamos la calidad de las imputaciones mediante la raíz del error cuadrático medio (para las observaciones categóricas no obstante carece de sentido esta medida). a = as.numeric(completos[is.na(data)]) b = as.numeric(datosoriginales[is.na(data)]) rmse = sqrt(mean((a-b)^2)) Otro ejemplo también con variables de distinto tipo muy sencillo es el siguiente: data(nhanes2) head(nhanes2) imputation = mice(nhanes2, m=10, maxit=100, method = rep("fastpmm", 4)) completos = complete(imputation, action = 5) print(completos)  Descubrimiento de temáticas Asignación latente de Dirichlet Formalización de la generación Inferencia y estimación usando muestreo de Gibbs colapsado Formalización del aprendizaje Ejemplos de aplicación Análisis de perfil en Twitter Temática según autores clásicos 8. Minería de texto 8.1 Descubrimiento de temáticas En este capítulo abordaremos un problema relativamente reciente de la Inteligencia Artificial. Supongamos que tenemos una gran cantidad de documentos que analizar, con el mismo o diferentes propósitos. Una de las problemáticas que nos puede surgir entre tanta información es la de buscar temáticas comunes entre diversos documentos. Esto puede tener mucha utilidad en diversos campos, pero especialmente en motores de búsqueda, donde las páginas (especialmente de noticias), quedan indexadas según una serie de palabras clave. Éstas serían posteriormente agrupadas automáticamente por el mismo motor para proporcionar información de una manera ordenada. En definitiva, nuestro problema explicado de manera informal es el siguiente: dados una serie de documentos de texto cualesquiera, determinar si existen k temáticas relativas a esos documentos. Este número k de temáticas viene fijado de antemano. En los últimos años se han desarrollado multitud de técnicas para intentar solventar este problema. Para este trabajo vamos a centrarnos en una llamada Latent Dirichlet Allocation (que traduciremos por Asignación latente de Dirichlet), cuya estimación de parámetros puede ser solucionada por técnicas MCMC. Nótese que estas técnicas no son las únicas para resolver el problema, pero como veremos se prestan a ser las más fáciles de implementar en este caso. La asignación latente de Dirichlet (la llamaremos LDA por comodidad a partir de ahora) está implementada en varios lenguajes de programación, siendo la más fácil de utilizar una de Python. Ésta, además implementa dicha técnica usando el muestreo de Gibbs (concretamente una versión llamada muestreo de Gibbs colapsado), por lo que resulta especialmente adecuada para ejemplificar. 8.2 Asignación latente de Dirichlet Antes de comenzar a explicar de manera formal en qué consiste dicha técnica, vamos a ver cómo funciona de manera constructiva: LDA supone que los documentos son mixturas de temáticas, y que las temáticas a su vez son mixturas de palabras. Es decir, supongamos en primer lugar que queremos generar un documento con respecto a estas suposiciones: 70 Capítulo 8. Minería de texto 1. Decidiríamos en primer lugar cuántas palabras N va a tener un documento, por ejemplo de acuerdo con una distribución de Poisson. 2. Decidimos una distribución para las temáticas, por ejemplo de acuerdo con una distribución Dirichlet. 3. Ahora generamos cada palabra del documento de la manera siguiente: a) Decidimos la temática de la palabra, escogiendo una de acuerdo a una distribución multinomial, conjugada con la Dirichlet anterior. b) Una vez escogida la temática, escogemos la palabra de acuerdo a la distribución multinomial de la temática. 8.2.1 Formalización de la generación Definamos lo siguiente: Una palabra es la unidad básica de dato discreto, definida como un miembro de vocabulario indexado entre 1,...,V . Representamos las palabras como vectores unitarios, que tienen un solo componente igual a 1 y el resto igual a 0. La v -ésima palabra del vocabulario queda representada por wv=1 solo para v. Un documento es una secuencia de Npalabras denotadas por w= (w1,...,wN). Un cuerpo es una colección de Mdocumentos denotados por D=w1,...,wM De manera formal, para generar un documento lo que haríamos sería lo siguiente: 1. Escoger N∼Poisson(ε) 2. Escogemos θ∼Dir(α) 3. Para cada una de las Npalabras wi: a) Escogemos una temática según zn∼Multinomial(θ) b) Escogemos una palabra wicon probabilidad P(wi|zi,β). Realizaremos algunas hipótesis de simplificación para continuar. Para empezar y como ya hemos comentado, la dimensionalidad k de la distribución Dirichlet queda fijada de antemano (y por tanto el número de temáticas a encontrar). En segundo lugar, las probabilidades de cada una de las palabras quedan parametrizadas por una matriz β de dimensiones k×V , donde βi j =P(wj=1|zj=1) . En un principio trataremos esta matriz como algo fijado que tendremos que estimar eventualmente. Por último, la suposición de que los documentos tienen un número de palabras de acuerdo con una distribución de Poisson puede ser ignorada completamente. 8.3 Inferencia y estimación usando muestreo de Gibbs colapsado Hasta ahora hemos aprendido cómo generar documentos de acuerdo con este modelo. Hemos empezado así porque LDA es un modelo generativo, en el sentido de que los documentos vienen generados según esa serie de hipótesis para luego realizar el aprendizaje de manera más sencilla. En primer lugar explicaremos de manera similar a como hemos hecho antes el procedimiento de aprendizaje de manera informal, y luego lo formalizaremos de manera más detallada. Nuestro caso ahora se centra en que tenemos una serie de documentos a los cuales queremos extraer temáticas: De manera aleatoria, para cada cada palabra en cada documento, asignamos una temática. Esta asignación ya nos da una representación de la distribución de las temáticas y de las palabras (si bien nada buenas). Para mejorarlas: • Para cada palabra en cada documento y para cada temática calcular 1) P(temticat|documento d) = la proporción de palabras en el documento d que están asignadas a la temática t y 2) P(palabra w|temtica t) = la proporción de asignaciones a la temática t debidas a esta palabra w. 8.3 Inferencia y estimación usando muestreo de Gibbs colapsado 71 • Reasignar a la palabra una nueva temática de acuerdo con P(temticat|documento d)× P(palabra w|temtica t)para todas las temáticas t Después de repetir este procedimiento un número elevado de veces, alcanzaremos un estado estable en el que las asignaciones apenas cambian. En este estado ya podemos usar las asignaciones para estimar las mixturas en cada documento (simplemente contando las palabras asignadas a cada temática en cada documento) y las palabras asociadas a cada temática (contando el total de palabras asignadas a cada temática). 8.3.1 Formalización del aprendizaje Dados D documentos conteniendo T temáticas expresadas bajo W palabras únicas, podemos representar P(w|z) como un conjunto de T distribuciones multinomiales φ sobre las W palabras, de tal manera que P(w|z=j) = φ(j) w y equivalentemente, representamos P(z) como D distribuciones multinomiales sobre las T temáticas, de tal manera que para un documento cualquiera P(z=j) = θ(d) j. Nuestra estrategia para descubrir las temáticas pasa por explorar la distribución a posteriori de las temáticas sobre las palabras P(z|w) . Una vez evaluada, podremos obtener estimaciones tanto de φ como de θ . (Estas dos últimas pueden ser representadas como matrices por conveniencia). Utilizamos el muestreo de Gibbs para muestrear de esta distribución a posteriori, para ello, definimos lo siguiente: wi|zi,φ(zi)∼Multinomial(φ(zi)) φ∼Dirichlet(β) zi|θ(di)∼Multinomial(θ(di)) θ∼Dirichlet(α) Donde α y β son hiperparámetros para controlar las distribuciones a priori de φ y θ . Estas distribuciones a priori son conjugadas con la Multinomial, permitiendo calcular la distribución conjunta P(w,z) = P(w|z)P(z) . Si en esta distribución integramos primero sobre φ obtendremos: P(w|z) = Γ(Wβ) Γ(β)WTT ∏ j=1 ∏wΓ(n(w) j+β) Γ(n(.) j+Wβ)(8.1) donde n(w) j es el número de veces que una palabra w ha sido asignada a la temática j en el vector de asignaciones z. Si en el segundo térmnino integramos sobre θ, obtendríamos: P(z) = Γ(Tα) Γ(α)TDD ∏ d=1 ∏jΓ(n(d) j+α) Γ(n(d) .+Tα)(8.2) donde n(d) j es el número de veces que una determinada palabra de un documento d ha sido asignada al documento j. Por tanto, nuestro objetivo sería evaluar la distribución a posteriori: P(z|w) = P(w,z) ∑zP(w,z)(8.3) Desgraciadamente, esta distribución no puede ser calculada directamente, ya que su cálculo implica la evaluación de una distribución de probabilidad sobre un espacio de estados discretos muy amplios. No obstante, podemos usar el muestreo de Gibss, para obtener la expresión, 72 Capítulo 8. Minería de texto mediante cancelación de términos de las dos distribuciones anteriores, obteniendo la siguiente distribución totalmente condicionada: P(zi=j|z−i,w)∝n(wi) −i,j+β n(.) −i,j+Wβ n(di) −i,j+α n(di) −i+Tα (8.4) donde n. −i es un conteo sin contar la asignación actual de zi . Este resultado es bastante intuitivo, ya que de la expresión anterior el primer término expresa la probabilidad de wi bajo la temática j , y el segundo ratio expresa la probabilidad de la temática j sobre el documento di . Una vez obtenida esta distribución totalmente condicionada, el algoritmo procede como sigue: Inicializamos todas las palabras de los documentos (realmente sólo aquellas incluídas entre las W del vocabulario) con valores entre 1 y T . A continuación asignamos temáticas a las palabras de acuerdo a las probabilidades obtenidas con la expresión anterior. Lo más adecuado, es que en cada iteración dichas probabilidades sean halladas únicamente con aquellas palabras cuya asignación haya sido realizado anteriormente. Una vez que la cadena ha corrido un número suficientemente alto de veces, llegará a un estado suficientemente estable (o estacionario). Con un conjunto de muestras de la distribución a posteriori ya tendríamos resuelto el problema de asignación. No obstante, si estuviéramos interesados en conocer la forma de los parámetros distribucionales podríamos estimarlos mediante dicha muestra de la siguiente manera: ˆ φ(w) j=n(w) j+β n(.) j+Wβ (8.5) ˆ θ(d) j=n(d) j+α n(d) .+Tα(8.6) Escogiendo los parámetros α,βy el número de temáticas T El algoritmo descrito en esta sección podría extenderse al caso en el que α y β fueran a su vez distribuciones dependientes de más hiperparámetros y muestrear de ellos, pero realizarlo realmente contribuye a un mayor coste computacional sin obtener mejoras de rendimiento. La elección de estos dos parámetros es importante a la hora de obtener distintos resultados: en particular, un β más grande produce un menos número de temáticas más amplias. Dados valores fijos a estos hiperparámetros se trataría de escoger un número de temáticas apropiado, lo que realmente es un problema de selección de modelos bayesiano. Lo natural en estadística bayesiana es evaluar la distribución a posteriori de los modelos dados los datos, y escoger aquel que de mayor densidad de probabilidad. En nuestro caso, nuestros datos son las palabras en los documentos, luego necesitaríamos evaluar la verosimilitud P(w|T). De nuevo, este es un problema intratable desde el punto de vista computacional, pero podríamos aproximar dicha distribución tomando la media harmónica de un conjunto de valores de P(w|z,T) cuando z se muestrea de la distribución a posteriori P(z|w,T) . Estas últimas pueden hallarse usando 8.1. Una vez aproximada, podemos utilizar esta cantidad para escoger el número de temáticas adecuadamente. 8.4 Ejemplos de aplicación En esta sección nos dedicaremos a aplicar este algoritmo en algunas situaciones interesantes en las que tenga sentido. Desde el punto de vista práctico estaremos utilizando el paquete lda disponible en PyPI , tanto bajo Python2 como Python3 . En este caso utilizaremos esta última implementación por evitar problemas con la autenticación via oauth2 . 8.4 Ejemplos de aplicación 73 8.4.1 Análisis de perfil en Twitter Con el auge de las nuevas redes sociales, un nuevo mundo de información textual se abre. Esta información es fácilmente accesible para cualquiera con una cuenta de por ejemplo Twitter (con algunas limitaciones) desde cualquier lenguaje de programación moderno. En nuestro caso, lo que vamos a realizar con el siguiente código es un análisis de perfil de una determinada cuenta de esta red social. Concretamente, estaríamos interesados en conocer el contenido en cuanto a temáticas de los tweets de una determinada persona. Centrémonos en analizar de qué habla Elon Musk, CEO de Tesla Motors. Para ello necesitamos en primer lugar una serie de claves de acceso a la API de Twitter, las cuales son obtenibles fácilmente a través de la página web. Comencemos el análisis. Ejercicio 8.1 #!/usr/bin/env python # -*- coding: utf-8 -*- """ Latent Dirichlet Allocation Twitter Crawler @author: JJimenez """ import numpy as np access_token = "xxxxx-xxxxxx" access_token_secret = "xxxxx" consumer_key = "xxxxxxxxxxxxxxxxxxxxxxxxx" consumer_secret = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" import tweepy auth = tweepy.OAuthHandler(consumer_key, consumer_secret) auth.set_access_token(access_token, access_token_secret) api = tweepy.API(auth) public_tweets = api.user_timeline('elonmusk', count = 1000) texts = [] for tweet in public_tweets: texts.append(tweet.text) from sklearn import feature_extraction paravocab = feature_extraction.text.CountVectorizer(stop_words = 'english', max_features = 50).fit(texts) vocab = paravocab.vocabulary_ vocab = vocab.keys() 80 BIBLIOGRAFÍA [MR14] Jean-Michel Marin y Christian P Robert. Bayesian Essentials with R. 2014, página 296. ISBN: 978-1-4614-8686-2. DOI: 10.1007/978-1-4614-8687-9 .URL: http://link.springer.com/10.1007/978-1-4614-8687-9 . [Mcm05] Data-driven Mcmc. “MCMC Tutorial at ICCV Motivation for MCMC”. En: October (2005). [MJG12] Morphometric Mcmc, Leif T Johnson y Charles J Geyer. “Morphometric MCMC (mcmc Package Ver. 0.9)”. En: (2012), páginas 1-20. [Mcm+05] Trans-dimensional Mcmc y col. “References • Peter Green , Reversible jump Markov chain Monte Carlo , in “ Highly Structured Stochastic Systems ”, 2003 • Paskin & Thrun , Robotic Mapping with Polygonal Outline Model Selection • Most common case : inference on p ( x | z ), x continuous p”. En: October October (2005), páginas 1-12. [Ml] An Ml. “Mixture models ( Ch . 16 ) Mixture models ( cont ’ d ) Bayesian estimation ( cont ’ d ) Mixture models - missing data”. En: 1 (). [Ric+14] Author Richard y col. “Package ‘ BayesFactor ’”. En: (2014). [Sah00] Sujit Sahu. “Markov Chain Monte Carlo ( MCMC ) Introduction”. En: August (2000). [VA13] Sunay Vaishnav y Primary Advisor. “A Markov Chain Monte Carlo based approach to Image Segmentation”. En: (2013). [WG00] Michael D Ward y Kristian Skrede Gleditsch. “Location , Location , Location : An MCMC Approach to Modeling Spatial Context with Categorical Variables in the Study and Prediction of War 1”. En: (2000). [Zhu+05] Song-chun Zhu y col. “Markov Chain Monte Carlo for Computer Vision (SLIDES)”. En: October October (2005). [ZS02] Zhuowen Tu y Song-Chun Zhu. “Image segmentation by data-driven markov chain monte carlo”. En: IEEE Transactions on Pattern Analysis and Machine Intelligence 24.5 (2002), páginas 657-673. ISSN: 01628828. DOI: 10.1109/34.1000239 .