Es dudoso que el género humano logre crear un enigma que el mismo ingenio humano no resuelva
Edgar Allan Poe
Como parte del curso de simulación que se imparte en el ITAM, se estudian diferentes metodologías para generar números seudo-aleatorios de diversas distribuciones de probabilidad. Entre los métodos que se revisan, se incluyen algunas de las versiones más elementales de la familia de métodos Monte Carlo basados en cadenas de Markov (MCMC) que son fundamentales para la estimación de integrales en el contexto Bayesiano.
Sin embargo, la amplitud de aplicaciones de los métodos MCMC puede ayudar a resolver problemas prácticos que en primera instancia no son evidentes en el estudio de estos métodos. El objetivo de esta nota es mostrar una aplicación que puede despertar mayor interés en estos métodos entre los estudiantes.
Introducción
El método de Monte Carlo basado en cadenas de Markov fue desarrollado por Nicholas Metropolis en el laboratorio de Los Álamos en 1953 como parte de la necesidad de resolver integrales multidimensionales en la simulación del comportamiento de un líquido y su fase gaseosa. Durante un largo tiempo, estos métodos fueron usados principalmente por químicos y físicos . Wilfred K. Hastings generalizó el algoritmo de Metropolis y el resultado de sus investigaciones fue publicado en su artículo de 1970 [5]. Hastings escribe: ’’When I learned how easy it was to generate samples from high dimensional distributions using Markov chains, I realised how important this was for Statistics, and I devoted all my time to this method and its variants which resulted in the 1970 paper’’.
Un caso especial del algoritmo de Metropolis-Hastings fue introducido por los hermanos Stuart Geman y Donald Geman en 1984, y se conoce como el ’’Gibss Sampler’’, nombrado así en honor a Josiah Williard Gibbs, un físico americano del siglo XIX que junto con Maxwell y Boltzmann, crearon la mecánica estadística. Fue hasta principios de los años noventas que comenzó a usarse de manera sistemática en la estadística, particularmente en la estadística Bayesiana, gracias al artículo de Smith y Gelfand [4]. Un recuento muy interesante del desarrollo del pensamiento y métodos Bayesianos se puede encontrar en el libro de Sharon Bertsch [2].
La idea básica de los métodos basados en MCMC es muy simple: queremos simular observaciones $X_1,X_2,\ldots,X_n$ de una densidad de probabilidad $f$. Entonces se construye una cadena de Markov $\{X_i\}$, fácil de simular y cuya distribución estacionaria $\pi$ corresponda a la distribución objetivo que nos interesa, es decir $\pi=f$. El teorema ergódico es la base de la aplicación del algoritmo MCMC.Algoritmo de Metropolis-Hastings
El algoritmo de Metropolis-Hastings (MH) simula una cadena de Markov cuya distribución estacionaria es $f$ de la siguiente manera:
- Se comienza con un valor inicial $X_0$ para la cadena.
- En el $n-$ésimo paso, dado un valor $X_n=x$ de la cadena, se simula un nuevo valor candidato $Y$, de una densidad propuesta $q(Y|X_n)$, es decir $Y\sim q(\cdot|X_n)$
- El valor propuesto $Y$ se podrá considerar o no como parte de la cadena, con una probabilidad $\alpha$ definida por $\alpha = \min\left\{1, \frac{f(y)q(x|y)}{f(x)q(y|x)} \right\}$. Entonces la cadena se actualiza con $X_{n+1} = y$ con probabilidad $\alpha$ o se queda en el valor actual $X_{n+1}=X_n$ con probabilidad $1-\alpha$.
A continuación veremos un ejemplo de la aplicación del algoritmo de Metropolis-Hastings, en un problema de descifrar un mensaje encriptado. El siguiente ejemplo se basa en el artículo de Persi Diaconis [3] y el desarrollo que hace Robert Dobrow [1].
Definición del modelo
Sea $A$ un conjunto que representa un alfabeto con $n$ letras y un espacio. Por ejemplo, en el inglés se consideran $n=26$ letras más el espacio son 27. Ahora consideremos un mensaje $M_k$ que consta de k letras en $A$. En el español podrían ser 28 o 29. Sea $\sigma: A \to A$ una permutación o función de codificación. Hay $n!$ permutaciones. Para el inglés como se definió arriba hay $27!\approx 10^{28}$ posibles funciones de codificación. \[ \sigma(c_1,c_2,\ldots,c_n) = (c_{\sigma(1)},c_{\sigma(2)},\ldots,c_{\sigma(n)}).\] Consideremos un mensaje $M_k$ conformado por $k$ letras en $A$. El problema a resolver es el siguiente: si $M_k$ es un mensaje cifrado, queremos encontrar la función $\sigma^*$ que decodifica el mensaje para hacerlo legible. Para descifrar el mensaje, se necesita medir de alguna manera la ‘distancia’ entre el mensaje cifrado y el mensaje real. El mensaje real no es conocido, por lo que se requiere una medición indirecta. Esta medición puede basarse en la verosimilitud de que ciertas letras aparezcan de manera conjunta. Para estimar esta verosimilitud, se puede usar una muestra del lenguaje común, que típicamente se conoce como corpus. Dado un corpus asociado al alfabeto $A$, podemos construir una función de evaluación $($como una especie de ‘distancia’$)$ de la siguiente manera: sea $O$ la matriz de $(n+1)\times (n+1)$ que registra la frecuencia de transiciones de una letra a otra en el corpus, formando pares de letras. Por ejemplo, si el corpus es: ‘abracadabra’, entonces $O_{ab}=O_{br}=O_{ra}=2$, $O_{ac}=O_{ca}=1$, etc. A la función $\sigma$ le podemos asignar la evaluación: \[ {\mbox eval}(\sigma) = \prod_{i=1}^k O_{\sigma(c_i),\sigma(c_{i+1})} \] La función de evaluación para el mensaje $M_k$ lo definimos como el producto sobre todas las frecuencias de los pares sucesivos de letras $(\sigma(c_i),\sigma(c_{i+1}))$ en el texto descifrado con la $\sigma$ particular. Esta evaluación será mayor cuando frecuencias de pares sucesivos en el mensaje corresponden con aquellos que aparecen más frecuentemente en el texto de referencia. Entonces, codificaciones con altas evaluaciones son las candidatas para decodificar. El objetivo es encontrar la $\sigma^*$ con valor de evaluación máxima. A continuación, se considera un corpus basado en la obra completa de Jane Austen, que se encuentra disponible en el Proyecto Gutenberg y de la que se calcula la matriz de transiciones $O$. Dobrow y Goldfeather hicieron el cálculo de la matriz, sería ideal contar con una versión similar en español.set.seed(123) #fija una semilla para reproducir el ejercicio.
# La siguiente matriz está basada en la obra de Jean Austen y cuenta las transiciones
# encontradas entre letras consecutivas (idioma inglés) Fuente: Dobrow.
archivo <- "https://raw.githubusercontent.com/jvega68/Simulacion/master/datos/AustenCount.txt"
mat <- read.table (archivo,header= F)
row.names(mat) <- colnames (mat) <- c (letters,"space")
mat[1:3,1:3] #transiciones de a->a, a->b a->c, b->a, y así sucesivamente
a b c
a 117 6669 7227
b 2385 208 0
c 8527 5 1843
logmat <- log(mat + 1) #transforma a logaritmos para facilitar operaciones.
Podemos definir una distribución de probabilidad que sea proporcional al score sobre el espacio de permutaciones, y el problema de Montecarlo consiste en muestrear de esta distribución:
\[\pi_{\sigma} = \frac{eval(\sigma)}{\sum_{\psi\in \mathcal{C}}eval(\psi)}\] Es importante notar que el denominador de esta densidad no se puede calcular, pues es la suma de $n!$ componentes. Pero para aplicar el algoritmo de Metropolis-Hastings sólo basta conocer los cocientes de la forma $\pi_{\sigma}/\pi_{\psi} = eval(\sigma)/eval(\psi)$, l cual simplificará mucho el procedimiento. evaluacion <-function(permutacion){
# calcula score del mensaje codificado con permutacion dada.
p <- 0
# Para cada par de letras en el mensaje decodificado
# busca la matriz de transición para la probabilidad de ese par
for (i in 1:(nchar(mensaje)-)){
p <- logmat[charIndex}(substr(permutacion, i, i)),
charIndex(substr(permutacion, i+1,i+1))]
}
# regresa la suma de estas probabilidades
return(p)
}
Elección del kernel de transición
Podemos ejecutar una caminata aleatoria en el conjunto de permutaciones de la siguiente manera: dada una permutación $\sigma$, la transición a una permutación propuesta $\sigma^*$ se da tomando dos letras al azar y cambiando los valores que $\sigma$ asigna a esas letras. Con este método de transposiciones aleatorias se contruye un kernel de transición simétrico, obteniendo como cociente de Hastings la siguiente expresión: \[ \alpha(\sigma,\sigma^*)= \frac{\pi_{\sigma^*} q(\sigma^*|\sigma)}{\pi_{\sigma}q(\sigma|\sigma^*)} = \frac{eval(\sigma^*)}{eval(\sigma)} \]Algoritmo
Para hacer las búsquedas de los pares en la matriz $O$ se indexan a través de sus códigos ASCII, para lo que se usan dos funciones: ascii y charIndex.ascii <- function(char){
# ascii(char) devuelve el numeral ascii de un caracter
strtoi (charToRaw(char),16L) #obtiene valor ascii crudo
}
charIndex <- function(char){
# esta función toma un caracter y regresa su 'char value'
# definido como a=1, b=2, ..., z=26, space=27 en la matriz
# esto es para que corresponda a la posición en la matriz diccionario
aValue <- ascii(char)
if (aValue == 32) {
27 #regresa 27 si es un espacio
} else {
aValue - 96 #ascii define "a" como 97, así que reescalamos restando 96}
}
Las evaluaciones ya hechas se van guardando en una lista de environments para no tener que recalcularlas cada vez.
Ejemplo
Consideremos el siguiente mensaje. Para simplificar, se han eliminado los signos de puntuación y la diferencia entre minúsculas y mayúsculas:
coincidences in general are great stumbling blocks in the way of that class of thinkers who have been educated to know nothing of the theory of probabilities that theory to which the most glorious objects of human research are indebted for the most glorious of illustrations edgar allen poe the murders in the rue morgue morpheus
#Este es el mensaje codificado. Hay que descodificarlo
mensaje <- "coincidences in general are great stumbling blocks in the way of that
class of thinkers who have been educated to know nothing of the theory of probabilities
that theory to which the most glorious objects of human research are indebted for the
most glorious of illustrations edgar allen poe the murders in the rue morgue morpheus"
decrypt <- function(perm, fa){
# Desencripta código de acuerdo al score actual fa
out <- perm
# para cada caracter en el mensaje, decodifica de acuerdo al score actual fa
for (i in 1:nchar(mensaje)){
CharInd <- charIndex(substr(perm,i,i))
if (charInd < 27) {
# cambia el i-ésimo caracter al determinado por fa
substr(out,i,i) <- rawToChar(as.raw(fa[charInd] + 96))
}
}
return(out)
}
codemess <- decrypt(mensaje,sample(1:26)) # codemess tiene el mensaje encriptado, con alguna permutación
codemess
[1] "nuxlnxcjlnjb xl kjljhot ohj khjog bgzistxlk stunwb xl gej fod ur geog ontobb ur gexlwjhb feu eoaj sjjl jcznogjc gu wluf lugexlk ur gej gejuhd ur mhusosxtxgxjb ogeog gejuhd gu fexne gej iubg ktuhxuzb usyjngb ur eziol hjbjohne ohj xlcjsgjc ruh gej oiubg ktuhxuzb ur xttzbghogxulb jckoh ottjl muj gej izhcjhb xl gej hzj iuhkzj iuhmejzb"
Entonces la simulación queda del siguiente modo:
[1] ’’300’’
[2] ’’wsouwodeuwea ou geuenir ine gneit atlycroug crswfa ou the miv sp thit wriaa
sp thoufena mhs hibe ceeu edlwited ts fusm usthoug sp the thesnv sp jnscicorotoea
thit thesnv ts mhowh the ysat grsnosla sckewta sp hlyiu neaeinwh ine oudected psn
the ysat grsnosla sp orrlatnitosua’’
[1] ’’600’’
[2] ’’cioucofeuces ou geuenal ane gneat strbploug plicms ou the way id that class
id thoumens whi have peeu efrcatef ti muiw uithoug id the theiny id jnipapolotoes
that theiny ti whoch the bist glinoirs ipkects id hrbau neseanch ane oufeptef din
the bist glinoirs id ollrstnatoius’’
[1] ’’900’’
[2] ’’cioucofeuces ou geuenal ane gneat strmploug plicbs ou the way id that class
id thoubens whi have peeu efrcatef ti buiw uithoug id the theiny id qnipapolotoes
that theiny ti whoch the mist glinoirs ipkects id hrmau neseanch ane oufeptef din
the mist glinoirs id ollrstnatoius’’
[1] ’’1200’’
[2] ’’coincifences in general are great stumpling plocbs in the way od that class
od thinbers who have peen efucatef to bnow nothing od the theory od qropapilities
that theory to which the most glorious opkects od human research are infeptef dor
the most glorious od illustrations’’
[1] ’’1500’’
[2] ’’coincidences in general are great stumpling plocbs in the way of that class
of thinbers who have peen educated to bnow nothing of the theory of jropapilities
that theory to which the most glorious opkects of human research are indepted for
the most glorious of illustrations’’
Podemos ver que después de un número de simulaciones el mensaje comienza a ser legible. Sin embargo, el número de iteraciones necesarias puede ser muy variable. Otro tema que requiere más estudio es entender el proceso de convergencia de las cadenas de Markov generadas y que puede ser motivo de otro artículo.
Conclusiones
En este artículo, se hizo una introducción muy general sobre los métodos MCMC, y se revisó un ejemplo de aplicación muy concreta en el ámbito de los modelos en criptografía.
Referencias
1Dobrow, Robert P. Introduction to Stochastic Processes with R. Wiley, 2016. 2 Bertsch McGrayne, SharonBibliografia
| 1. | ↑ | Dobrow, Robert P. Introduction to Stochastic Processes with R. Wiley, 2016. |
| 2. | ↑ | Bertsch McGrayne, Sharon
|
| 3. | ↑ | Diaconis, Persi.
``The Markov chain Monte Carlo revolution’’,
|
| 4. | ↑ | Gelfand, Alan E. y Smith, Adrian F. M.
``Sampling-Based Approaches to Calculating Marginal Densities’’,
|
| 5. | ↑ | Hastings, W. K.
``Monte Carlo sampling methods using Markov chains and their applications’’,
|