Estadística II (Grado en Matemáticas)

Añado a la página de cursos del blog la mayor parte del material que he estado utilizando los dos últimos cursos en la asignatura Estadística II, optativa en el Grado en Matemáticas de la UAM. Dado que es un material de clase que tal vez requeriría más tiempo para ser revisado es posible que aún queden errores e imprecisiones por lo que debe utilizarse con precaución.

Se trata de una segunda asignatura de estadística después de otra obligatoria, en la que se estudia la teoría básica de inferencia. El núcleo central de esta segunda asignatura es el modelo de regresión lineal. Previamente se dan dos temas, uno sobre la distribución normal multivariante que proporciona la base para poder demostrar muchos de los resultados del curso, y otro sobre contrastes de hipótesis no paramétricas (bondad de ajuste, homogeneidad e independencia). El curso termina con una introducción a los métodos de clasificación supervisada en la que se tratan dos clasificadores lineales: el basado en el modelo de regresión logística y la regla de Fisher.

Una variable aleatoria que casi siempre toma valores menores que su media

Si la distribución de una variable aleatoria es asimétrica la media puede dar una impresión equivocada de los valores que toma la variable. Dada una probabilidad   tan cercana a uno como se quiera, en esta entrada veremos como dar un ejemplo sencillo de variable aleatoria tal que , es decir, si elegimos un valor de cercano a uno, toma valores que son casi siempre menores que su valor esperado . Posteriormente veremos que este ejemplo da lugar a un estimador insesgado que se comporta de una manera inesperada. He encontrado el ejemplo en  este blog.

El ejemplo

Sea una variable normal estándar y sea . Definimos , es decir, es una v.a. con distribución log-normal de parámetros y . La esperanza de esta variable es . Por lo tanto, dado arbitrario,

,

para un valor de suficientemente grande.

Estimación de la media de una variable log-normal

Supongamos que es una muestra de v.a.i.i.d. con distribución log-normal de parámetros y y que queremos estimar la media de la población. Consideramos los dos estimadores alternativos siguientes:

• La media muestral: este es un estimador insesgado que, además, no requiere ningún conocimiento de cuál es la distribución que siguen los datos.
• Sabiendo que la población es log-normal, se puede usar también el  estimador plug-in que resulta de la expresión , donde . Esto lleva a , donde es la varianza muestral de . Este estimador tendrá un pequeño sesgo que no debería ser muy importante para tamaños muestrales grandes.

¿Cómo funcionan estos estimadores? Para responder a esta pregunta he simulado muestras muy grandes (de tamaño ) para distintos valores de . Para cada muestra he calculado el cociente entre cada uno de los dos estimadores anteriores y el valor del parámetro. Los resultados son los que muestra el gráfico:

Los puntos negros son realizaciones de un estimador insesgado, mientras que los blancos corresponden a un estimador que sí tiene sesgo. Así que los resultados parecen ser justo al revés de lo que dice la intuición. El estimador insesgado prácticamente siempre toma valores menores que el parámetro. Sin embargo, su sesgo es cero porque en ocasiones con probabilidad muy pequeña toma un valor muy grande que compensa las numerosas veces en las que el parámetro se infraestima. En la simulación ha ocurrido una vez (el punto negro de la parte superior derecha).

A pesar de no ser insesgado, el estimador plug-in es mucho mejor que la media muestral, para estimar la media de la población. Esto no sorprende porque incorpora el conocimiento de cuál es la distribución de los datos. La media muestral en este caso funciona mal porque hereda el comportamiento de la variable, que está casi siempre por debajo de su media.

Este es el código de R usado en la simulación:

set.seed(100)
n <- 10^6
sigma <- seq(0.1, 5, 0.1)

theta1 = numeric(length(sigma))
theta2 = numeric(length(sigma))
theta = numeric(length(sigma))

for (i in 1:length(sigma)){
z <- rnorm(n)
x <- exp(sigma[i]*z)
theta[i] = exp(sigma[i]^2/2)
theta1[i] = mean(x)
theta2[i] = exp(sd(log(x))^2/2)
}

plot(sigma, theta1/theta, pch=16, xlab='sigma',
ylab='estimador / parámetro', xlim=range(sigma))
abline(h=1)
points(sigma, theta2/theta)
legend('topleft', c('media muestral','plug in'), pch=c(16,1))

Dos formas de abordar los problemas de probabilidad

La lectura del problema planteado en este tuit y las correspondientes respuestas que obtuvo me llevó a pensar en qué significa realmente resolver un problema de probabilidad y qué consecuencias para nuestra docencia puede tener la respuesta a esta pregunta.

Today’s probability question: pic.twitter.com/z6osB0x7y7

— Michael Lopez (@StatsbyLopez) 14 de abril de 2017

Algunas repuestas hicieron un sencillo cálculo de probabilidad para hallar la respuesta:

@StatsbyLopez Difference between two uniform distributions is a triangular distribution.
Percentage of area in peak sixth is (5/6 + 6/6)/6 = 11/36 pic.twitter.com/kr5uQqYoih

— Darren O’Shaughnessy (@RankingSw) 14 de abril de 2017

Sin embargo, el problema se puede resolver (con el grado de precisión que se quiera) mediante una simulación que requiere solo una línea en R:

@StatsbyLopez 30.55%

mean(abs(runif(1e7, 0, 2) – runif(1e7, 0, 2)) < 1 / 3)

— David Robinson (@drob) 14 de abril de 2017

A efectos prácticos ambos enfoques son igualmente válidos. Con la solución analítica se consigue una comprensión más profunda de algunos aspectos (por ejemplo, se sabe cómo es la distribución de la diferencia de dos distribuciones uniformes). Sin embargo,  la respuesta basada en simulación es más potente que la analítica en varios sentidos. Si cambiamos la distribución uniforme por otra más complicada tendríamos que rehacer completamente los cálculos, pero bastaría cambiar una palabra en la línea de código. Además, si no sabemos si el tiempo en el que David y Emily llegan al Rec Center sigue una distribución uniforme, pero hemos registrado el tiempo al que han llegado durante 10 días (y suponemos que los tiempos de llegada son independientes), podríamos muestrear de la distribución empírica lo que conecta la solución del problema de probabilidad con una metodología estadística importante (el bootstrap).

# Generamos los datos
set.seed(2) # para reproducir resultados
n = 10
david = runif(n, 0, 2)
emily = runif(n, 0, 2)

#Muestreamos de las dist. empíricas
mean(abs(sample(david, 1e7, rep=TRUE) - sample(emily, 1e7, rep=TRUE)) < 1/3)

 

Para los datos generados por el código anterior, el resultado es 0.30 aprox., no lejos del valor obtenido en el caso en que conocemos que las variables tienen distribución uniforme.

Mi conclusión es que los métodos basados en simulación son muy importantes y creo que en el futuro lo serán aún más. Me parece que deberían ganar peso en los programas de nuestras asignaturas. Algunos autores han argumentado que la enseñanza de la inferencia mediante simulaciones debería constituir la base de los cursos habituales de introducción a la estadística [véase Tintle et al. (2015) y las muchas referencias que se citan en este artículo]. Copio del abstract de Tintle et al. (2015) las ventajas pedagógicas de esta metodología:

Termino con enlaces a tres entradas del blog relacionadas con esta. Las dos primeras plantean problemas complicados de resolver analíticamente, pero que resultan fáciles de abordar mediante simulación. La tercera corresponde a un artículo que escribí sobre las interacciones entre los métodos de simulación y la inferencia estadística:

• El problema de los cuatro puntos
• La distancia entre dos puntos aleatorios en un círculo
• Simulación e inferencia estadística

 

Referencia

Tintle, N., Chance, B., Cobb, G., Roy, S., Swanson, T., y VanderStoep, J. (2015). Combating Anti-Statistical Thinking Using Simulation-Based Methods Throughout the Undergraduate Curriculum. The American Statistician, 69, 362-370.

El problema de los cuatro puntos

En esta entrada se plantea un problema de probabilidad complicado de resolver analíticamente pero cuya solución puede aproximarse fácilmente por métodos de Montecarlo (simulación).

En 1864, Sylvester propuso la primera versión de un problema de geometría estocástica que se conoce con el nombre de problema de los cuatro puntos:

Al seleccionar cuatro puntos aleatoriamente en el plano, demostrar que vale 1/4 la probabilidad de que el conjunto convexo más pequeño que los contiene (su cierre convexo) sea un triángulo.

Un ejemplo de situación en la que ocurre el suceso cuya probabilidad queremos calcular es el siguiente:Por el contrario, si el punto negro no perteneciera al triángulo definido por los tres puntos rojos no habría ocurrido el suceso en el que estamos interesados ya que el cierre convexo sería un cuadrilátero.

Tal y como lo planteó Sylvester el problema es ambiguo mientras no se dé un significado preciso a la frase seleccionar cuatro puntos aleatoriamente en el plano. Esto explica que diferentes autores obtuvieran diferentes probabilidades, y diferentes a 1/4, en función del significado que atribuían a esa frase. Esto es lo que ocurre también en la más conocida paradoja de Bertrand.

Posteriormente se planteó el problema de forma completamente precisa: sea un cuerpo convexo (es decir, un conjunto convexo y compacto con interior no vacío). Si se seleccionan cuatro puntos aleatoria y uniformemente sobre , ¿cuál es la probabilidad de que el cierre convexo de los cuatro puntos sea un triángulo?

La solución depende de y se puede demostrar que  para cualquier , con lo que la probabilidad siempre está entre 0.29 y 0.33 más o menos.

Una solución aproximada mediante simulación  

A continuación vamos a aproximar la solución del problema para cuatro puntos seleccionados en el cuadrado unidad mediante una pequeña simulación. La siguiente función de R genera cuatro puntos en el cuadrado unidad y devuelve TRUE o FALSE en función de si el cierre convexo es un triángulo o no, y este experimento se lleva a cabo B veces.

sylvester.cuadrado = function(B){
replicate(B, length(chull(cbind(runif(4), runif(4)))) == 3)
}

La obtención del cierre convexo con R mediante el comando chull  (así como la aplicación del cierre convexo en ecología) la hemos tratado en una entrada anterior del blog.

Hemos usado la función anterior para seleccionar (10000 veces) 4 puntos en el cuadrado unidad. Posteriormente calculamos la proporción de veces en las que el cierre convexo resultó ser un triángulo:

# Ejemplo
set.seed(100) # para reproducir los resultados
B = 10000
resultado = sylvester.cuadrado(B)
sum(resultado) / B

En el ejemplo anterior, a partir de 10000 realizaciones obtenemos la aproximación . En 1867, Woolhouse probó que el valor exacto para es . Puede comprobarse que si es una transformación afín no singular, entonces . Por lo tanto, la probabilidad para cualquier rectángulo es también 11/35.

Ejercicio: escribir una función de R similar a la anterior para aproximar el valor de cuando es un círculo (o elipse) y cuando es un triángulo. Comprobar que, aparentemente, la probabilidad para un círculo es menor que para un cuadrado, y ésta es menor que para un triángulo.

El problema variacional 

Algo más tarde, en 1885, Sylvester se interesó por el problema de qué forma debería tener para que alcanzara su valor máximo y su valor mínimo. Conjeturó que el valor mínimo se alcanza para el círculo y el máximo para el triángulo. Una prueba rigurosa de esta conjetura no se obtuvo hasta 1917 y se debe a Blaschke.

Los resultados mencionados en esta nota han sido tomados de Pfiefer (1989). Allí se puede encontrar que los valores exactos de para el triángulo y el círculo son, respectivamente, 1/3 y , valores que permitirán comprobar si la solución del ejercicio propuesto más arriba es correcta y que, junto con la solución del problema variacional, justifican las cotas para que mencionamos más arriba. También se pueden encontrar las principales ideas involucradas en las demostraciones de los resultados.

Referencia

Pfiefer, R.E. (1989). The historical development of JJ Sylvester’s four point problem. Mathematics Magazine, 309–317.

¿Son incorreladas dos variables aleatorias independientes cuya covarianza es cero?

Sean e dos variables aleatorias independientes cuya covarianza es igual a cero. ¿Podemos decir que el coeficiente de correlación lineal entre ellas también vale cero?

En esta entrada se presenta un ejemplo que demuestra que, estrictamente hablando, la respuesta a la pregunta es negativa. La razón es que para que el coeficiente de correlación esté bien definido hay que imponer condiciones a la existencia de momentos de las variables. Además, estas condiciones son ligeramente más restrictivas que las necesarias para que exista la covarianza.

Se trata de una cuestión sobre todo técnica, pero es necesario tenerla en cuenta al aplicar el coeficiente de correlación a variables con colas pesadas como las que aparecen, por ejemplo, en datos económicos. Y en cualquier caso hay que tenerla en cuenta para hacer afirmaciones matemáticamente correctas.

El ejemplo

Como punto de partida consideramos v.a.i.i.d. con distribución normal estándar. Definimos e . Obviamente e son v.a. independientes. Además  e verifican las siguientes propiedades:

1. .
2. , ya que es la inversa de una v.a. con distribución chi-cuadrado con 4 grados de libertad. (Se aplica la fórmula que aparece en esta entrada de wikipedia.)
3. Dado que tiene distribución con 1 y 4 grados de libertad  también se verifica (de nuevo, se puede consultar la fórmula en wikipedia) y, por lo tanto, .
4. Como consecuencia de los tres puntos anteriores .
5. Sin embargo, ya que la varianza de la inversa de una chi-cuadrado con grados de libertad solo es finita si . Por lo tanto el coeficiente de correlación entre e no está definido.

He tomado el ejemplo anterior de Mukhopadhyay (2010). En el artículo se presenta también un ejemplo aún más sencillo en el que las variables son independientes pero ni siquiera la covarianza está definida. Basta tomar e , con lo que tiene distribución de Cauchy y no es finita.

Referencias

Mukhopadhyay, N. (2010). When finiteness matters: counterexamples to notions of covariance, correlation, and independence. The American Statistician, 64, 231-233.

 

VI jornada estadística en la UAM

El próximo día 24 de febrero, viernes, tendremos en el Departamento de Matemáticas de la UAM una nueva jornada estadística, con diferentes charlas de miembros del departamento y de algunos profesores invitados de otras universidades. Este es el programa completo de la jornada:

Una (dulce) receta con R

El brazo de gitano (llamado en los países anglosajones Swiss roll) es un conocido bizcocho relleno y enrollado en forma de cilindro. En esta entrada se describe cómo generar una muestra aleatoria en que tiene la forma de este dulce, y cómo representar estos datos en un gráfico tridimensional.

Este ejercicio de simulación no es una simple curiosidad, puesto que estos datos o similares se utilizan en la literatura estadística sobre manifold learning y métodos no lineales de reducción de la dimensión para ilustrar las propiedades de los métodos. Un ejemplo es el artículo original acerca de Isomap [Tenenbaum et al. (2000)], que es un conocido método de reducción no lineal de la dimensión.

Según Izenman (2008), página 617, las ecuaciones de la variedad con forma de brazo donde toman valores las observaciones son las siguientes: , y donde y . En nuestra receta se generan   realizaciones de una variable aleatoria uniforme en el intervalo y otras realizaciones (independientes de las anteriores)  de una variable aleatoria con distribución uniforme en el intervalo . Insertando los valores obtenidos en las ecuaciones se obtienen los puntos.

Para representar los datos en un gráfico tridimensional dinámico he usado el comando plot3d de la librería rgb. El código en R y el gráfico tridimensional dinámico (los puntos se pueden girar y hacer más grandes o pequeños con el ratón) se pueden ver aquí. El resultado de nuestra receta es el siguiente:

Imaginemos que queremos identificar dos clusters diferentes en una muestra de datos en esta variedad. Los métodos clásicos de clustering ( medias, por ejemplo) se basan en distancias entre las observaciones que no tienen en cuenta la estructura no lineal de la variedad en la que viven, por lo que en general no van a dar buenos resultados. Adoptar un punto de vista adecuado es crucial. Los gráficos siguientes corresponden a los mismos datos contemplados desde dos lugares diferentes:

Mientras que en el gráfico de la izquierda no se distingue apenas la existencia de grupos bien diferenciados en los datos, en el de la derecha aparecen dos clusters de forma totalmente natural. Resulta por lo tanto necesario disponer de métodos que permitan detectar la estructura no lineal de la variedad en la que viven los datos. Este problema puede ser bastante complicado si la dimensión del espacio ambiente es mucho mayor que la de la variedad a la que pertenecen realmente las observaciones.

Volviendo al auténtico brazo de gitano dulce con el que empezábamos la entrada, en esta página se cuentan varias historias que explican su curioso nombre en España. Según wikipedia en otros países de habla hispana se le llama de otras formas, también curiosas: brazo de reina (Chile y Colombia), niño envuelto (México), arrollado (Argentina) o pionono (Perú).

Referencias

Izenman, A. J. (2008). Modern multivariate statistical techniques. Springer.

Tenenbaum, J. B., De Silva, V. y Langford, J. C. (2000). A global geometric framework for nonlinear dimensionality reduction. Science, 290, 2319-2323.