Una introducción al paquete CARET

El paquete CARET (classification and regression training, Kuhn (2016)) incluye una serie de funciones que actúan como interfaz para decenas de métodos complejos de clasificación y regresión. Utilizar este paquete en lugar de las funciones originales de los métodos presenta dos ventajas:

  • Permite utilizar un código unificado para aplicar reglas de clasificación muy distintas, implementadas en diferentes paquetes.
  • Es más fácil poner en práctica algunos procedimientos usuales en problemas de clasificación. Por ejemplo, hay funciones específicas para dividir la muestra en datos de entrenamiento y datos de test o para ajustar parámetros mediante validación cruzada.

He tenido que usar este paquete para elaborar algunos ejemplos de un artículo y, de paso, he escrito una introducción muy básica con un resumen de lo que he aprendido. La podéis encontrar en este enlace.

Una “chuleta” que resulta útil para recordar lo básico del uso de Caret se puede encontrar en este enlace.

Anuncios
Publicado en estadística, R | Etiquetado | Deja un comentario

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

La Facultad de Ciencias de la Universidad Autónoma de Madrid (UAM) y la Unidad de Bioestadística del Instituto IMDEA Alimentación han organizado conjuntamente el curso de formación continua en Estadística Aplicada con R (tercera edición) que se va a impartir entre septiembre y diciembre de 2017 en Madrid de forma presencial

El programa está compuesto por seis módulos que se pueden cursar de forma independiente:

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

Todos los módulos se impartirán en las aulas de informática de la Facultad de Ciencias de la UAM.  Más información en el siguiente enlace: http://tinyurl.com/ydgo79bf .

Publicado en estadística, R | Etiquetado | Deja un comentario

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.

Publicado en estadística | Etiquetado , , , , | Deja un comentario

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  p\in (0,1) tan cercana a uno como se quiera, en esta entrada veremos como dar un ejemplo sencillo de variable aleatoria X tal que \mathbb{P}(X < \mathbb{E}(X)) > p, es decir, si elegimos un valor de p cercano a uno, X toma valores que son casi siempre menores que su valor esperado \mathbb{E}(X). 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 Z una variable normal estándar y sea \sigma>0. Definimos X=\exp(\sigma Z), es decir, X es una v.a. con distribución log-normal de parámetros 0 y \sigma. La esperanza de esta variable es \mathbb{E}(X)=\exp(\sigma^2/2). Por lo tanto, dado p\in (0,1) arbitrario,

\mathbb{P}(X < \mathbb{E}(X)) = \mathbb{P}(\sigma Z < \sigma^2/2) = \mathbb{P}(Z < \sigma/2) > p,

para un valor de \sigma suficientemente grande.

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

Supongamos que X_1,\ldots,X_n es una muestra de v.a.i.i.d. con distribución log-normal de parámetros 0 y \sigma y que queremos estimar la media \mu = \exp(\sigma^2/2) 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 \mu = \exp(\sigma^2/2), donde \sigma^2 = \mbox{Var}(\log(X)). Esto lleva a \hat{\mu} =\exp(S^2/2), donde S^2 es la varianza muestral de \log X_1,\ldots,\log X_n. 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 10^6) para distintos valores de \sigma. 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:

LogNormal

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))
Publicado en estadística, probabilidad | Etiquetado | Deja un comentario

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.

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

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:

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:

Abstract-Tintle et al-2015-TAS

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:

 

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.

Publicado en estadística, probabilidad | Etiquetado , | Deja un comentario

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:ejemplosylvesterPor 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 S\subset \mathbb{R}^2 un cuerpo convexo (es decir, un conjunto convexo y compacto con interior no vacío). Si se seleccionan cuatro puntos aleatoria y uniformemente sobre S, ¿cuál es la probabilidad p(S) de que el cierre convexo de los cuatro puntos sea un triángulo?

La solución depende de S y se puede demostrar que \frac{35}{12\pi^2} \leq p(S)\leq \frac{1}{3} para cualquier S, 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 S=[0,1]^2 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 p(S)\approx 0.302. En 1867, Woolhouse probó que el valor exacto para S=[0,1]^2 es p(S)=11/35\approx 0.3056. Puede comprobarse que si T es una transformación afín no singular, entonces p(T(S))=p(S). 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 p(S) cuando S es un círculo (o elipse) y cuando S 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 S para que p(S) 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 p(S) para el triángulo y el círculo son, respectivamente, 1/3 y 35/(12\pi^2), 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 p(S) 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.

Publicado en probabilidad, R | Etiquetado , , | Deja un comentario

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

Sean X e Y 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 Z_1,\ldots,Z_5 v.a.i.i.d. con distribución normal estándar. Definimos X=Z_1^2 e Y=(Z^2_2+Z^2_3+Z^2_4+Z^2_5)^{-1}. Obviamente X e Y son v.a. independientes. Además X e Y verifican las siguientes propiedades:

  1. E(X)=E(Z_1^2)=1.
  2. E(Y)=1/2, ya que Y 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 4XY tiene distribución F con 1 y 4 grados de libertad  también se verifica E(4XY)=2 (de nuevo, se puede consultar la fórmula en wikipedia) y, por lo tanto, E(XY)=1/2.
  4. Como consecuencia de los tres puntos anteriores Cov(X,Y)=E(XY)-E(X)E(Y)=0.
  5. Sin embargo, Var(Y)=\infty ya que la varianza de la inversa de una chi-cuadrado con n grados de libertad solo es finita si n>4. Por lo tanto el coeficiente de correlación entre X e Y 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 X=Z_1 e Y=Z_2^{-1}, con lo que XY=Z_1/Z_2 tiene distribución de Cauchy y E(XY) no es finita.

Referencias

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

 

Publicado en estadística, probabilidad | Etiquetado , | Deja un comentario