VIII Jornada Estadística UAM

VIIIJornadaUAM

El próximo día 29 de marzo de 2019 a las 10:00 tendremos en la sala 520 del Departamento de Matemáticas (módulo 17 de la Facultad de Ciencias) la VIII 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:

Captura

Anuncios
Publicado en estadística | Etiquetado , | 2 comentarios

Acerca de jamovi

¿Qué es jamovi?

Jamovi es un programa gratuito para realizar cálculos estadísticos. Su apariencia es la de una hoja de cálculo, pero es una interfaz gráfica de R (en inglés, una GUI, Graphical User Interface). Cada vez que elegimos una opción del menú de análisis, en realidad se ejecuta código de R para llevar a cabo los cálculos, sin embargo R está encapsulado dentro del programa y el usuario no se entera de ello, si así lo desea.

jamovi.JPG

Diseñado para competir con SPSS

El programa parece estar diseñado para competir con SPSS. El aspecto general de los cuadros de diálogo es muy similar y no tiene problema para abrir los ficheros de datos .sav de SPSS como si fueran propios.

jamovi1

Creo que tiene algunos rasgos generales que mejoran a SPSS:

  • Cuando se elige alguna opción se realiza el cálculo de forma inmediata sin necesidad de confirmacion. Esta característica lo hace más ágil que SPSS.
  • Los resultados van apareciendo en una ventana y al pedir que se salve el análisis se genera un fichero con extensión .omv en el que se guarda todo junto: los datos, los resultados y las opciones que han dado lugar a estos resultados. Esto es mejor que lo que hace SPSS (que guarda los resultados en un fichero .spo diferente y no guarda las opciones) ya que permite reproducir los cálculos de manera mucho más fiable. Por otra parte, usando el botón derecho del ratón se guardan los resultados numéricos parciales (en .pdf o .html) y los gráficos en diversos formatos (.pdf, .png, etc.)
  • Si se modifican los datos, se vuelven a rehacer los cálculos automáticamente. Esto también es más ágil que en SPSS.

Posible software para un curso de introducción a la estadística

El menú tiene bastantes menos funciones que SPSS, pero suficientes para un curso de introducción a la estadística. En cualquier caso las funciones se pueden ampliar a través de una serie de módulos que se instalan fácilmente. La idea es que progresivamente los usuarios puedan añadir más módulos, a la manera de los paquetes de R.

La versión que yo he manejado incluye, entre otras funciones:

  • Descripción de datos numérica y gráfica básica. El diagrama de dispersión requiere instalar el módulo scatr (hubiera sido mejor incluirlo directamente en el menú Descriptives). Con unos pocos clics se puede obtener, por ejemplo, este diagrama de dispersión con datos de dos grupos que incluye las dos rectas de mínimos cuadrados y diagramas de cajas para las distribuciones marginales:

dispersion

  • Contrastes e intervalos básicos: t de Student, test \chi^2 de bondad de ajuste… En los contrastes es fácil elegir si se desea un llevar a cabo un contraste bilateral o unilateral (opción que SPSS tampoco tiene).
  • Diseño de experimentos y análisis de la varianza incluyendo algún procedimiento no paramétrico.
  • Regresión lineal y logística, incluyendo gráficos de residuos, técnicas de diagnóstico, etc.
  • Componentes principales y análisis factorial.
  • Análisis de supervivencia, añadiendo el módulo Death Watch (de nombre un tanto tétrico).

Mi impresión es que se trata de un proyecto con posibilidades aunque aún no se sabe cómo va a evolucionar y si al final la comunidad de usuarios será grande. Lo he estado usando (alternándolo con SPSS) en un curso de introducción a la estadística que acabo de impartir y la experiencia ha sido positiva: el manejo básico se aprende en 10 minutos, lo que se aprende es útil para trabajar después con SPSS si es necesario y los estudiantes lo pueden instalar sin problemas de licencia en sus ordenadores.

Opciones que he echado de menos

He echado de menos algunas opciones:

  • No es posible aún editar los gráficos, nos tenemos que conformar con las opciones elegidas por los programadores.
  • Se puede instalar un editor de R, de manera que podamos ejecutar código de R simple. Esto está muy bien si queremos que los alumnos se vayan familiarizando con R, pero el editor es muy rudimentario y no he podido guardar el código en un fichero .R.
  • Algunas técnicas básicas en clasificación no supervisada o supervisada: análisis de clústeres (k medias por ejemplo) o la regla de clasificación de Fisher.

Recursos

Algunos enlaces útiles para el que se anime a probar:

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

Bendición de la dimensionalidad y clasificación perfecta de datos funcionales

Cuando la dimensión de los datos es muy alta (se miden para cada individuo muchas variables) muchas regiones del espacio muestral están vacías de datos por lo que algunos métodos estadísticos fallan. Un ejemplo son los métodos que requieren optimizar funciones complicadas que dependen de muchas variables, ya que el espacio a explorar buscando el óptimo es demasiado extenso. Este fenómeno se conoce con el nombre de maldición de la dimensionalidad.

Sin embargo medir más variables también significa tener más información, lo que para la mayoría de los propósitos debería ser positivo. En esta entrada se da un ejemplo de bendición de la dimensionalidad en un problema de clasificación de datos funcionales. Es la explicación que subyace al fenómeno de “clasificación perfecta” descrito por Delaigle y Hall en 2012. A su vez, como veremos, este fenómeno tiene conexiones con la teoria de espacios de Hilbert con núcleo reproductor (RKHS).

Clasificación supervisada

Se dispone de dos muestras de funciones procedentes de dos posibles poblaciones. El problema consiste en clasificar una nueva función, independiente de las muestras anteriores y cuya procedencia se desconoce, en una de las dos posibles poblaciones. Cuestiones relacionadas con este problema ya las hemos tratado en este blog anteriormente.

Si en lugar de funciones se dispone de vectores de p variables el problema es muy conocido. Para simplificar, supongamos que la población 1 es normal p-dimensional de media m y matriz de covarianzas \Sigma mientras que la población 0 es normal p-dimensional de media 0 y la misma matriz de covarianzas. Además, suponemos que a priori la probabilidad de las poblaciones es la misma. En este caso la regla óptima de clasificación es la regla de Fisher: la observación x se clasifica en la población 1 si w'x-w'm/2>0, donde w=\Sigma^{-1}m. Además, la mínima probabilidad de error posible (el llamado error Bayes) es 1-\Phi(\Delta/2), donde \Phi es la función de distribución normal estándar y \Delta = (m'\Sigma^{-1}m)^{1/2} es la distancia de Mahalanobis entre los centros de las dos poblaciones. Esta expresión del error es bastante intuitiva: si \Delta=0 las poblaciones son iguales y la mejor probabilidad de error posible es 1/2 (la que se obtendría tirando una moneda). A medida que \Delta aumenta, las dos poblaciones se van separando y el error Bayes va disminuyendo. La mínima probabilidad de error es función decreciente de \Delta.

El efecto de considerar cada vez más variables: dos posibilidades

Al añadir más variables la distancia de Mahalanobis \Delta aumenta y por lo tanto el error Bayes disminuye. En este sentido tenemos una bendición de la dimensionalidad. ¿Cuándo al añadir un número suficiente de variables se obtiene una clasificación casi perfecta (el error Bayes converge a cero) y cuándo no (error Bayes por encima de cierto nivel estrictamente positivo)? Ese es esencialmente el cálculo que Delaigle y Hall llevaron a cabo en 2012.

Supongamos que la matriz de covarianzas tiene rango p y llamamos \lambda_1>\lambda_2>\cdots>\lambda_p>0 a sus autovalores (que suponemos todos diferentes) y v_1,\ldots,v_p a los correspondientes autovectores ortonormales. Podemos diagonalizar la matriz \Sigma de la siguiente forma \Sigma = \sum_{i=1}^p \lambda_i v_iv'_i. También podemos escribir el vector de medias en función de la base v_1,\ldots,v_p, m=\sum_{i=1}^p m_i v_i. Como consecuencia, la distancia de Mahalanobis también se puede expresar como \Delta^2 = m'\Sigma^{-1} m =\sum_{i=1}^p m_i^2/\lambda_i. Si añadimos una nueva variable tal que la suma anterior aumenta mucho, entonces se producirá una gran descenso de la probabilidad de error de clasificación.

Si añadimos una sucesión infinita de variables, entonces alcanzaremos una clasificación perfecta si y solo si la serie \sum_{i=1}^\infty m_i^2/\lambda_i es divergente (es decir cuando la distancia de Mahalanobis crece a infinito, \Delta\to\infty). Si \sum_{i=1}^\infty m_i^2/\lambda_i<\infty, por muchas variables que añadamos siempre habrá una probabilidad de error de clasificación por encima de cierto valor estrictamente positivo y no podrá haber clasificación perfecta.

Aunque puede parecer un poco artificial eso de “añadir una sucesión infinita de variables”, cuando los datos son funcionales esa es exactamente la situación. Si suponemos que el espacio en el que viven las funciones es por ejemplo L_2[0,1] podemos desarrollar cada función en términos de la base ortonormal formada por las autofunciones del operador de covarianzas. De esta forma, cada función (y también la función media) se puede identificar con la sucesión de coeficientes de la función en el desarrollo de la serie. Por ejemplo, la función de medias es m=\sum_{i=1}^\infty m_i v_i, donde m_i =\langle m, v_i\rangle. Así pues, se produce clasificación perfecta de datos funcionales si y solo si la relación entre la media y la estructura de covarianzas es tal que \sum_{i=1}^\infty m_i^2/\lambda_i diverge. Esta es la conclusión a la que llegaron Delaigle y Hall. Con un número finito de variables es imposible alcanzar clasificación perfecta salvo en casos degenerados (que algún autovalor se anule, lo que ocurre con probabilidad cero en la práctica).

Si no hay clasificación perfecta, la regla óptima en el caso de datos funcionales es la extensión inmediata de la fórmula de Fisher, es decir, clasificar x en la población 1 si \langle w, x\rangle - \langle w, m\rangle /2 > 0, donde w = \sum_{i=1}^\infty (m_i/\lambda_i) v_i. Esta es la regla del centroide en la terminología de Delaigle y Hall. Para que este clasificador tenga sentido hay que asumir w\in L_2[0,1], es decir, \sum_{i=1}^\infty m^2_i/\lambda^2_i < \infty. Esta condición es estrictamente más fuerte que suponer que no hay clasificación perfecta (\sum_{i=1}^\infty m^2_i/\lambda_i < \infty).

La conexión con los espacios de Hilbert con núcleo reproductor

En la primera parte de un artículo publicado este año [Berrendero, Cuevas y Torrecilla (2018)] hemos reinterpretado estas expresiones en términos de la teoría RKHS mejorando así los resultados de Delaigle y Hall. Es posible demostrar que hay clasificación perfecta de datos funcionales si y solo si la función de medias m no pertenece al RKHS definido por la función de covarianzas del proceso cuyas trayectorias queremos clasificar. Si la función de medias sí pertenece al RKHS y no hay clasificación perfecta en el artículo damos una expresión completamente general de la regla óptima que no se puede expresar en términos “elementales”, es decir, no se puede expresar en términos de productos escalares en L_2[0,1]. De hecho, en la regla óptima pueden intervenir integrales estocásticas, como ocurre en el caso del movimiento browniano que vamos a ver un poco más adelante.

La conexión con la dicotomía de Feldman-Hájek

Dos medidas de probabilidad son equivalentes si dan probabilidad cero a exactamente los mismos sucesos y son mutuamente singulares si existe un suceso que tiene probabilidad cero para una medida y probabilidad uno para la otra. La llamada dicotomía de Feldman-Hájek afirma que dos medidas gaussianas sobre espacios muy generales, y en particular en L_2[0,1], son siempre o equivalentes o mutuamente singulares.

En nuestro artículo observamos que la condición de clasificación perfecta de datos funcionales que hemos visto es equivalente a que las medidas de probabilidad correspondientes a las dos poblaciones entre las que queremos clasificar sean mutuamente singulares. Si queremos clasificar una función en este caso y existe un conjunto de funciones con probabilidad uno para una población y probabilidad cero para la otra, basta identificar si la función está en este conjunto o no para clasificarla sin probabilidad de error. En la práctica puede ser difícil llevar a cabo esta identificación.

Un ejemplo: el movimiento browniano

Consideramos dos muestras de movimientos brownianos estándar, en la primera de ellas -la de color azul en el gráfico- la media es cero (consiste en realizaciones independientes de B(t),\, t\in[0,1]) mientras que en la segunda -la de color rojo- hay una tendencia determinista m (realizaciones independientes de m(t) + B(t),\, t\in [0,1]). El siguiente gráfico muestra un ejemplo de esta situación. Dada una nueva curva cuyo color es desconocido, tenemos que clasificarla como procedente de la población de las rojas o de las azules:

ejemplo-fda1

La función de covarianzas del browniano es K(s,t)=\min\{s,t\} y se puede demostrar que el RKHS correspondiente es el conjunto de funciones m(t) con m(0)=0 y tales que existe m'\in L_2[0,1] tales que m(t)=\int_0^t m'(u)du. Si la tendencia está en este conjunto no hay clasificación perfecta. Esta es una condición de suavidad. Por ejemplo, si la tendencia es una función lineal, entonces no hay clasificación perfecta.

Puede demostrarse que la regla óptima en este ejemplo consiste en clasificar una función X(t) en la población roja si \int_0^1 m'(u)dX(u) - (1/2)\int_0^1m'(u)^2du >0. En la expresión de la regla óptima aparece la integral estocástica de Itô.
En el gráfico anterior la tendencia es lineal positiva m(t)=2t. En este caso la regla óptima es extremadamente simple puesto que \int_0^1 m'(u)dX(u) - (1/2)\int_0^1m'(u)^2du=2X(1)-2, se clasifica una función X(t) en la población con tendencia si x(1)>1, la trayectoria acaba por encima de la línea horizontal del gráfico siguiente:

Bayes-Brown

Es llamativo en este ejemplo que lo único que necesitamos saber de una trayectoria para clasificarla de forma optima es el nivel en que acaba. De hecho, en el artículo damos una condición suficiente sobre la función media m para que la regla de clasificación óptima solo dependa de un numero finito de puntos de la trayectoria. Este resultado se toma posteriormente como base para un procedimiento de selección de variables en clasificación de datos funcionales.

Para simplificar la discusión anterior he supuesto que todos los parámetros de las poblaciones son conocidos. En la práctica al añadir más variables -y en particular si tratamos con datos funcionales- los problemas de estimación que hay que resolver son cada vez más complejos, y es en este punto donde la dimensionalidad comienza a mostrar su lado maldito.

Referencias

Berrendero, J.R., Cuevas, A. y Torrecilla, J.L. (2018). On the use of reproducing kernel Hilbert spaces in functional classification. Journal of the American Statistical Association, 113, 1210-1218. Preprint | Material suplementario

Delaigle, A., y Hall, P. (2012). Achieving near perfect classification for functional data. Journal of the Royal Statistical Society: Series B, 74, 267-286.

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

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

Si X es una variable aleatoria con función de distribución continua F, entonces la variable U=F(X) tiene distribución uniforme en el intervalo [0,1]. Por ejemplo, si tenemos observaciones independientes  x_1,\ldots, x_n procedentes de una distribución normal estándar y llamamos \Phi a la función de distribución \mbox{N}(0,1), entonces \Phi(x_1),\ldots, \Phi(x_n) es una muestra de observaciones independientes de una uniforme en [0,1].

¿Existe alguna transformación similar para convertir vectores aleatorios de dimensión p (con una distribución continua cualquiera) en vectores uniformes en el (hiper)cubo [0,1]^p? 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 (X_1,X_2) con distribución continua. La transformación de Rosenblatt consiste en calcular (u_1,u_2) = T(x_1, x_2), donde

  • u_1 = F_1(x_1), donde F_1(x_1) es la función de distribución marginal de X_1.
  • u_2 = F_2(x_2|x_1), donde F_2(x_2|x_1) es la distribución de X_2 condicionada a X_1=x_1.

Se puede demostrar que el vector (U_1,U_2) = T(X_1,X_2) tiene distribución uniforme en el cuadrado unidad [0,1]^2.

Por ejemplo, consideremos el caso en que (X_1,X_2) tiene distribución normal bidimensional con vector de medias \mu=(\mu_1,\mu_2) y matriz de covarianzas

\Sigma=\left( \begin{array}{ll} \sigma_1^2 & \sigma_{12} \\ \sigma_{21} & \sigma_{2}^2 \end{array} \right).

Por las propiedades de la distribución normal X_1 tiene distribución normal de media \mu_1 y varianza \sigma_1^2 mientras que la distribución de X_2 condicionada a X_1=x_1 es también normal con media \mu_{2.1} = \mu_2 + \sigma_{12}(x_1-\mu_1)/\sigma_1^2 (la recta de regresión) y varianza \sigma_{2.1}^2 = \sigma_2^2 - \sigma_{12}^2/\sigma_1^2. De acuerdo con estas fórmulas, la transformación que tenemos que hacer es:

  • u_1 = \Phi[(x_1-\mu_1)/\sigma_1], donde \Phi es la función de distribución normal estándar.
  • u_2 = \Phi[(x_2-\mu_{2.1})/\sigma_{2.1}].

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 F_0. 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:

EjemploMuestras

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:

NoNormal

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 p la transformación es la siguiente:

  • u_1 = F_1(x_1), donde F_1(x_1) es la función de distribución marginal de X_1.
  • u_2 = F_2(x_2|x_1), donde F_2(x_2|x_1) es la distribución de X_2 condicionada a X_1=x_1.
  • u_p = F_p(x_p|x_1,\ldots,x_{p-1}), donde F_p(x_p|x_1,\ldots,x_{p-1}) es la distribución de X_p condicionada a X_1=x_1,\ldots,X_{p-1}=x_{p-1}.

Una última observación es que en realidad hay p! 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 Statistics23, 470-472.

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

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 n pasos y el otro m 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?

DosBorrachos

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 n pasos de longitudes X_1,X_2,\ldots,X_n 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 D_n a la distancia entre el punto final y el origen tras n pasos, entonces el teorema de Rayleigh  proporciona la siguiente sencilla fórmula para n>2:

P(D_n < 1) = \frac{1}{n+1}

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 m y n pasos con n+m>2 y llamamos D_m y D_n a las correspondientes distancias entre los puntos finales e iniciales de cada una de ellas, entonces:

P(D_n < D_m) = \frac{m}{n+m}

El teorema de Rayleigh es un caso particular para m=1 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):

 

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():

trayectoria

Referencia

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

 

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

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:

optima.PNG

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:

optima-prod

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:

optima2.PNG

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 x en una de las poblaciones si y solo si f(x)>0, donde f(x)=x'Ax+2a'x+c, entonces para calcular los productos de hiperplanos basta diagonalizar la matriz B=aa'-cA, que depende de manera muy sencilla de la matriz A, el vector a y el escalar c . 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 B 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 B 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 ClassificationPreprint | doi

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

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:

GrandSlam

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:

shanghai

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:

biometrika

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.

 

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