Cómo transformar un vector aleatorio para que tenga distribución uniforme

Si es una variable aleatoria con función de distribución continua , entonces la variable tiene distribución uniforme en el intervalo . Por ejemplo, si tenemos observaciones independientes  procedentes de una distribución normal estándar y llamamos a la función de distribución , entonces es una muestra de observaciones independientes de una uniforme en .

¿Existe alguna transformación similar para convertir vectores aleatorios de dimensión (con una distribución continua cualquiera) en vectores uniformes en el (hiper)cubo ? En 1952 Murray Rosenblatt publicó un artículo de tres páginas en Annals of Mathematical Statistics en el que describe una transformación con esta propiedad.

En el apartado siguiente describimos con detalle la transformación en el caso bidimensional, en especial en el importante caso en que los datos tienen distribución normal bidimensional. Luego comentamos una posible aplicación.

La transformación de Rosenblatt en el caso bidimensional

Consideremos vectores bidimensionales con distribución continua. La transformación de Rosenblatt consiste en calcular , donde

• , donde es la función de distribución marginal de .
• , donde es la distribución de condicionada a .

Se puede demostrar que el vector tiene distribución uniforme en el cuadrado unidad .

Por ejemplo, consideremos el caso en que tiene distribución normal bidimensional con vector de medias y matriz de covarianzas

.

Por las propiedades de la distribución normal tiene distribución normal de media y varianza mientras que la distribución de condicionada a es también normal con media (la recta de regresión) y varianza De acuerdo con estas fórmulas, la transformación que tenemos que hacer es:

• , donde es la función de distribución normal estándar.
• .

Aplicación a los contrastes de bondad de ajuste

La transformación en el caso unidimensional se usa, por ejemplo, para demostrar que el estadístico del contraste de Kolmogorov-Smirnov tiene siempre la misma distribución bajo la hipótesis nula (es de distribución libre) cuando la distribución de los datos bajo la nula es continua. En este blog también hemos usado esta propiedad en una entrada anterior para resolver un problema suponiendo datos uniformes sin perder generalidad.

En el caso multivariante la propiedad también se puede usar para construir un contraste de bondad de ajuste a una distribución continua dada . Si aplicamos la transformación de Rosenblatt, el problema se reduce a saber contrastar si nuestras observaciones transformadas tienen distribución uniforme en un cubo.

El siguiente código de R es una implementación de esta transformación para datos normales bidimensionales:

library(MASS)
# Rosenblatt (p = 2)

rosenblatt = function(x, mu, sigma){
# x es una matriz de n x 2 datos
# mu es el vector de medias teorico de x
# sigma es la matriz de covarianzas teorica de x
mu2.1 = mu[2] + sigma[1,2] * (x[,1] - mu[1]) / sigma[1,1]
sigma2.1 = sigma[2,2] - sigma[1,2]^2/sigma[1,1]
u1 = pnorm((x[,1] - mu[1])/sqrt(sigma[1,1]))
u2 = pnorm((x[,2] - mu2.1)/sqrt(sigma2.1))
return(cbind(u1, u2))
}

#-------------------------------------------------------------

n = 1000
mu = c(2, 3)
sigma = matrix(c(1, 0.8, 0.8, 1), 2)
x = mvrnorm(n, mu, sigma)
colnames(x) = c("x1", "x2")
u = rosenblatt(x, mu, sigma)

layout(matrix(1:2, 1))
plot(x, pty="s", pch=16, cex=0.5)
plot(u, pty="s", pch=16, cex=0.5)

La segunda parte del código anterior genera una muestra aleatoria de datos normales bidimensionales con correlación 0,8 y la transforma en una muestra uniforme en el cuadrado unidad:

Si por ejemplo generamos datos normales, luego los elevamos a 1,2 (con lo que dejan de ser normales) y aplicamos la misma transformación anterior obtenemos:

Si detectamos la falta de uniformidad en el gráfico de la derecha (que parece clara por los puntos que se acumulan en la parte superior derecha del gráfico) esto permite rechazar la hipótesis de normalidad en los datos de la izquierda. Como hemos dicho antes, cualquier contraste de bondad de ajuste a una distribución continua bidimensional se reduce a saber contrastar la uniformidad en un cuadrado.

El código del último ejemplo es el siguiente:

n = 1000
mu = c(2, 3)
sigma = matrix(c(1, 0.8, 0.8, 1), 2)
x = mvrnorm(n, mu, sigma)^(1.2)
colnames(x) <- c('x1', 'x2')
u = rosenblatt(x, mu, sigma)

layout(matrix(1:2, 1))
plot(x, pty='s', pch=16, cex=0.5)
plot(u, pty='s', pch=16, cex=0.5)

 

El caso general

En el caso de dimensión arbitraria la transformación es la siguiente:

• , donde es la función de distribución marginal de .
• , donde es la distribución de condicionada a .
• …
• , donde es la distribución de condicionada a .

Una última observación es que en realidad hay transformaciones, ya que cada permutación de las variables da lugar a una transformación diferente. En el caso bidimensional podríamos haber comenzado transformando la segunda coordenada y luego la primera coordenada condicionando al valor de la segunda. Esto hubiera dado lugar a otra muestra transformada diferente, pero también uniforme en el cuadrado unidad.

Referencias

Rosenblatt, M. (1952). Remarks on a multivariate transformation. The Annals of Mathematical Statistics, 23, 470-472.

Anuncios

Historia de dos borrachos

Dos borrachos salen de un bar y dan pasos de un metro de longitud, cada paso en una dirección aleatoria e independiente de la dirección de los pasos anteriores. Si uno de ellos da pasos y el otro pasos, ¿cuál es la probabilidad de que el primero acabe más cerca del bar que el segundo? ¿Cuál es la probabilidad de que cualquiera de los dos acabe de nuevo en el bar (a distancia menor de un metro)?

A continuación se muestran dos trayectorias de 1000 pasos (en rojo) y 3000 pasos (en azul). El bar es el punto negro y los puntos finales de los dos paseos se han marcado en rojo y azul, respectivamente. En este caso el borracho que da más pasos ha acabado más cerca. ¿Cuál es la probabilidad de que esto ocurra?

Caminos o paseos aleatorios

El término paseo aleatorio (random walk) fue utilizado por primera vez por Karl Pearson en una carta a la revista Nature publicada en 1905 cuyo título era The problem of the random walk. Ese mismo año Einstein publicó también su trabajo sobre el movimiento browniano, que puede entenderse como límite de una sucesión de caminos aleatorios. El símil del borracho ya aparece en la propia carta de Pearson que, comentando un resultado de Rayleigh, acaba con la frase

 The lesson of Lord Rayleigh’s solution is that in open country the most probable place to find a drunken man who is at all capable of keeping on his feet is somewhere near his starting point!

El teorema de Rayleigh

Este resultado se puede usar para responder a la segunda de las preguntas que hacíamos al principio de la entrada. Consideramos caminos aleatorios que comienzan en el origen del plano y se mueven pasos de longitudes en direcciones aleatorias independientes con distribución uniforme. ¿Cuál es la probabilidad de que una trayectoria acabe en un lugar que dista menos de una unidad del punto de partida? Si llamamos a la distancia entre el punto final y el origen tras pasos, entonces el teorema de Rayleigh  proporciona la siguiente sencilla fórmula para :

Como consecuencia de este teorema la probabilidad de que dando 1000 pasos, el borracho acabe a menos de un metro del bar es 1/1001.

En Bernardi (2013) se proporciona una breve demostración de un resultado más general: si las longitudes de los pasos son variables aleatorias i.i.d. no negativas, consideramos dos trayectorias independientes de y pasos con y llamamos y a las correspondientes distancias entre los puntos finales e iniciales de cada una de ellas, entonces:

El teorema de Rayleigh es un caso particular para y pasos de longitud constante igual a 1.

En nuestro ejemplo inicial la probabilidad de que el borracho que da 3000 pasos acabe más cerca del bar que el que da 1000 es de 1/4. Las apuestas a favor de quedar más cerca del bar dando 1000 pasos que dando 3000 (odds ratio) están 3 a 1, lo cual es bastante intuitivo. Esta probabilidad no depende de la distribución de la variable aleatoria que representa la longitud de los pasos, siempre que sea la misma para los dos borrachos.

Código en R

El código en R utilizado para obtener las trayectorias de los caminos aleatorios anteriores es prácticamente el mismo que el que aparece en el siguiente tuit de Antonio S. Chinchón (@aschinchon) cuya lectura me ha sugerido esta entrada (las únicas diferencias son que he sustituido geom_polygon() por geom_path(), he elegido las direcciones con distribución uniforme continua en lugar de discreta y he añadido los puntos final e inicial de cada trayectoria):

 

library(ggplot2)

df <- data.frame(x=0, y=0)

for (i in 2:2500)
{
a <- sample(0:15,1)
df[i,1] <- df[i-1,1]+(cos(pi/8*a))
df[i,2] <- df[i-1,2]+(sin(pi/8*a))
}

ggplot(df, aes(x,y)) + geom_polygon() + theme_void()#rstats pic.twitter.com/L3gefSpWZf

— Antonio S. Chinchón (@aschinchon) August 16, 2018

En realidad esta entrada sirve como excusa para señalar lo elegantes y misteriosas que resultan las trayectorias cuando se sustituye geom_path() por geom_polygon(). Por ejemplo, esta es una trayectoria correspondiente a 2500 pasos cuando se usa geom_polygon():

Referencia

Bernardi, O. (2013). A Short Proof of Rayleigh’s Theorem with Extensions. The American Mathematical Monthly, 120, 362-364. Preprint.

 

Las componentes lineales de una regla de clasificación cuadrática

¿Cuándo una regla lineal es equivalente a una cuadrática?

Las reglas de clasificación cuadráticas son muy flexibles y permiten distinguir poblaciones que difieren no solo en las medias sino también en la estructura de covarianzas. De hecho, si queremos clasificar entre dos poblaciones  normales con matrices de covarianzas distintas, la regla óptima (o regla Bayes) es cuadrática.

Pero incluso en situaciones en las que la regla de clasificación óptima es cuadrática, es posible que no se pierda nada (en términos de error de clasificación) al sustituirla por una regla lineal más simple. Consideremos el ejemplo siguiente:

En este ejemplo tenemos dos muestras de entrenamiento procedentes de poblaciones normales con matrices de covarianzas diferentes. La regla de clasificación óptima es cuadrática y corresponde a la línea discontinua del gráfico. Está claro que si la cambiamos por una regla lineal no perdemos nada en error de clasificación y que la complejidad que introducen los términos cuadráticos es superflua en este caso.  Podemos usar una regla lineal, mucho más simple y fácil de interpretar.

Aproximación de reglas cuadráticas mediante productos de hiperplanos

En un artículo que se publicará proximamente hemos considerado el problema de cómo aproximar una regla cuadrática mediante productos de hiperplanos, con la idea de identificar situaciones como la del ejemplo anterior. La ventaja es que las aproximaciones solo dependen de combinaciones lineales de las variables, lo que permite definir reglas de clasificación más fáciles de interpretar. Además, es posible saber cuándo el producto es parecido a la función cuadrática original. En el ejemplo anterior el resultado de aplicar nuestro método a los datos disponibles es el siguiente:

En este caso solo uno de los dos hiperplanos del producto aproximante (las líneas continuas) es necesario para obtener una regla lineal con propiedades equivalentes a la regla cuadrática original. En otros casos pueden ser necesarios los dos hiperplanos o incluso otros productos de hiperplanos complementarios. Consideremos el ejemplo siguiente:

El gráfico de arriba a la derecha corresponde a la aproximación a la regla optima basada en un único producto de hiperplanos. Su error de clasificación va a ser bastante mejor que el obtenido mediante una regla lineal (por ejemplo Fisher) al precio de considerar dos funciones lineales en lugar de una. Incluso nos podemos acercar más a la regla cuadrática óptima si consideramos un segundo producto de hiperplanos (gráfico de abajo a la derecha). Para ello tendríamos que tener en cuenta cuatro combinaciones lineales de las variables originales.

La idea es que combinando productos de hiperplanos se obtiene una regla de clasificación lineal a trozos que combina las ventajas de las reglas cuadráticas (error de clasificación) y de las lineales (facilidad de interpretación). Los hiperplanos que intervienen en los productos constituyen la estructura lineal de la regla cuadrática.

El método se aplica a cualquier regla cuadrática, no necesariamente a las óptimas. Usando técnicas de regularización se han propuesto muchas reglas cuadráticas incluso en casos de alta dimensión en los que no es posible estimar las matrices de covarianzas. Un artículo reciente en el que se mencionan muchas de ellas es Qin (2018).

¿Es muy difícil determinar la estructura lineal de una regla cuadrática?

Si la regla cuadrática consiste en clasificar en una de las poblaciones si y solo si , donde , entonces para calcular los productos de hiperplanos basta diagonalizar la matriz , que depende de manera muy sencilla de la matriz , el vector y el escalar . Y esto para cualquier dimensión. De manera similar a como ocurre en el análisis de componentes principales, si muchos de los autovalores de son pequeños es posible aproximar bien la regla cuadrática con un número pequeño de productos de hiperplanos. De esta manera también se puede considerar nuestro método como una forma de reducir la dimensión del problema.

Si el lector tiene curiosidad por saber por qué la matriz es justo la que hace falta diagonalizar para encontrar los productos de hiperplanos aproximantes tendrá que leer nuestro artículo 🙂

Referencia

Berrendero, J.R. and Cárcamo, J. (To appear). Linear components of quadratic classifiers. Advances in Data Analysis and Classification. Preprint | doi

El Grand Slam de las revistas estadísticas

Según wikipedia el término Grand Slam procede de juegos de cartas de la familia del whist y en este contexto ya se usaba en el siglo XVII. Desde 1930 se usa en golf para referirse al hecho de ganar los cuatro principales torneos en la misma temporada. Al parecer el periodista Alan Gould  fue el primero en usar este término referido al tenis en un artículo de 1933 en el periódico The Reading Eagle:

Pero, ¿cuál sería el Grand Slam de las revistas de estadística? Si a los investigadores en estadística les preguntaran cuáles son las revistas más importantes, ¿cuál sería su respuesta?

Una encuesta

Recientemente se han publicado los resultados de 2018 del ranking de Shanghái. Junto con la lista de mejores universidades (ya suficientemente comentada en diversos foros) se publican también los resultados de una encuesta en la que se pregunta a científicos de distintas áreas, entre otras cosas, por las revistas de su especialidad que consideran más prestigiosas. Me ha parecido interesante mirar los resultados porque reflejan la importancia de las revistas de una forma más cualitativa que los controvertidos índices de impacto. En el área de estadística son cuatro las revistas que aparecen, lo que podríamos llamar el Grand Slam de las revistas estadísticas. Estos son los resultados:

Si comparamos con el ranking de revistas en función del índice de impacto hay ciertas discrepancias. Annals, JASA y Royal B son (todavía) Q1, aunque no son ciertamente las tres revistas con mayor índice de impacto. Por otra parte, hay muchos años en los que Biometrika no aparece como Q1, con lo que ello implica en procesos de acreditación, solicitudes de sexenios, etc.

Biometrika

Y sin embargo, a quién no le gustaría publicar en una revista fundada por Pearson, Weldon y Galton a principios del siglo XX y en la que han aparecido contribuciones fundamentales en la historia de la metodología estadística. Así cuenta David Cox los comienzos de Biometrika en un artículo publicado para celebrar su centenario:

En 1908 Gosset publicó en Biometrika su artículo sobre la distribución t de Student. Más de un siglo después miles de científicos siguen usando a diario esta distribución. Me cuesta imaginar un artículo que contribuya más al impacto de una revista científica.

 

¿Con qué probabilidad existe el estimador de máxima verosimilitud en regresión logística?

Un nuevo resultado sobre el estimador de máxima verosimilitud en regresión logística

En una entrada anterior de este blog se comentaba un conocido resultado sobre la existencia del estimador de máxima verosimilitud (EMV) en el modelo de regresión logística. El EMV existe si y solo si las nubes de puntos de las variables regresoras correspondientes a los dos grupos de la variable respuesta ( e ) no se pueden separar linealmente mediante un hiperplano  (la separación no tiene por qué ser estricta).

En un artículo publicado recientemente en arxiv, Emmanuel Candes y Pragya Sur han profundizado en esta situación cuando las variables regresoras tienen distribución normal. Si es el tamaño muestral y es el número de variables regresoras, en el artículo se investiga con qué probabilidad podemos esperar que exista un hiperplano separador (y por lo tanto el EMV no exista) cuando y , con  . Esta probabilidad es relevante en la práctica porque, dada la gran cantidad de datos a nuestra disposición, muchas veces se puede presentar la posibilidad de ajustar el modelo de regresión logística con un número elevado de variables regresoras.

El resultado es sorprendente ya que existe un brusco cambio de régimen en función de si está por encima o por debajo de cierto umbral. Si está por debajo, el EMV existe en el límite con probabilidad uno; mientras que si está por encima, el EMV existe con probabilidad cero. El umbral solo depende del valor del término independiente del modelo y del límite de las varianzas , donde es el vector de coeficientes de la regresión y es el vector de variables regresoras.

Una pequeña simulación

Tenía curiosidad por comprobar la relevancia empírica de este resultado y, en particular, si se notaba claramente el cambio de régimen para tamaños muestrales moderados. Para ello, he hecho una pequeña simulación y comparto aquí los resultados.

He considerado el caso más sencillo en el que  , es decir, la variable respuesta es independiente de las variables regresoras. En este caso, puede comprobarse que el umbral que determina el cambio de régimen es . Por lo tanto, si el número de variables regresoras está por debajo de la mitad del tamaño muestral debemos esperar probabilidades altas de existencia del EMV mientras que si está por encima se espera una probabilidad de existencia baja. Hay que decir que el umbral crítico es una función decreciente de y por lo tanto si fuese el umbral sería aún menor.

Para he generado muestras con variables regresoras independientes con distribución normal estándar, para y, en cada caso, he comprobado si son linealmente separables (doy más detalles al final sobre esta comprobación). Este experimento lo he repetido 100 veces. En el siguiente gráfico se muestra la proporción de veces en las que el EMV existe en función del número de variables regresoras:

Los resultados coinciden totalmente con lo esperado. Hasta una dimensión cercana a 100 el EMV existe siempre. En un entorno de 100 la proporción de veces en las que existe desciende bruscamente a cero y se mantiene en cero a partir de valores de la dimensión ligeramente por encima de 100.

A continuación he repetido el experimento para tamaños muestrales aún menores ( y ) con los resultados siguientes:

 

Vemos que ya para el efecto del cambio de régimen se nota claramente. Incluso para tamaños muestrales pequeños no es recomendable ajustar un modelo de regresión logística con un número de variables regresoras por encima del umbral que describen Candes y Pragya en su artículo. Desgraciadamente este umbral no es conocido porque depende en general de los valores de los parámetros que queremos estimar.

¿Qué ocurre para variables regresoras con otras distribuciones?

El teorema se cumple bajo condiciones más generales, al menos en el caso en que es independiente de . De hecho, bajo la condición hay resultados generales ya en los años sesenta debidos a Tom Cover [Cover, 1965]. A continuación aparece el resultado que he obtenido simulando variables regresoras independientes con distribución exponencial de media uno:

 

Comprobación de la separabilidad lineal

Para comprobar si existe un hiperplano separador he seguido el enfoque basado en optimización lineal propuesto por Konis (2007) Para ello, he usado la función separator incluida en el paquete de R safeBinaryRegression. Como la función está un poco escondida, yo la he copiado directamente en mi código, que se puede encontrar aquí por si alguien quiere repetir o modificar lo que yo he hecho.

Referencias

Cover, T. M. (1965). Geometrical and statistical properties of systems of linear inequalities with applications in pattern recognition. IEEE transactions on electronic computers, 3, 326-334.

Candes, E. y Pragya, S. (2018). The phase transition for the existence of the maximum likelihood estimate in high-dimensional logistic regression. https://arxiv.org/abs/1804.09753

Konis, K., B. (2007). Linear programming algorithms for detecting separated data in binary logistic regression models. PhD thesis, University of Oxford.

Curso de estadística aplicada con R (cuarta edición)

Los Departamentos de Matemáticas y Ecología de la Facultad de Ciencias de la Universidad Autónoma de Madrid organizan, en colaboración con el departamento de Bioestadística de GEICAM, la cuarta edición del curso Estadística aplicada con el software R. El programa está compuesto por diferentes módulos que pueden ser cursados en su totalidad o separadamente, dependiendo de los diferentes intereses y conocimientos de los alumnos. Los módulos se impartirán entre septiembre y noviembre de 2018 en la Facultad de Ciencias de la UAM:

1. Introducción a R (este módulo se impartirá en dos turnos, el 6-7 y el 20-21 de septiembre)
2. Métodos de regresión y análisis multivariante con R (3 a 5 de octubre)
3. Métodos de Regresión Avanzados para la Investigación en Ciencias Naturales con R (16 a 18 de octubre).
4. Estadística Aplicada a la Investigación Biomédica con R (22 a 26 de octubre).
5. Modelos Mixtos / Jerárquicos / Multinivel con R (20 a 22 de noviembre).
6. Técnicas Estadísticas de Data Mining con R (26 a 30 de noviembre).

Más información en este enlace.

VII Jornada Estadística UAM

El próximo día 23 de febrero de 2018, viernes, tendremos en la sala 520 del Departamento de Matemáticas (módulo 17 de la Facultad de Ciencias) la VII jornada estadística UAM, con diferentes charlas de miembros del departamento y de algunos profesores invitados de otras universidades. Este es el programa completo de la jornada: