Si considero un vector aleatorio \(\bf X\) con distribución uniforme, ¿la distancia de Mahalanobis de un punto a la esperanza de la distribución será cero porque la densidad uniforme es constante?.
No. En realidad, la distancia de Mahalanobis sólo utiliza la información de varianzas (dispersión) y correlaciones (relación lineal) del vector aleatorio, no toda la información contenida en la densidad. Así que en realidad \(d_M\) sólo toma en cuenta una geometría muy básica de los conjuntos de nivel de la densidad. El hecho de que \(d_M({\bf x},{\bf 0})=d_E({\bf y},{\bf 0})\), donde \(d_E\) es la distancia euclídea e \(\bf y\) es la estandarización multivariante de \(\bf x\), quiere decir que primero centramos \(\bf x\), luego lo medimos en unidades de desviación típica y, luego le quitamos su estructura de correlaciones y finalmente medimos la distancia euclídea del vector resultante, \({\bf y}\), a 0.
Cuando la densidad de \(\bf X\) es normal, entonces eso sí quiere decir que \(d_M({\bf x},{\boldsymbol\mu})\) es más pequeña para vectores \({\bf x}\) en los que la densidad \(f({\bf x})\) es más parecida a \(f({\boldsymbol\mu})\). Pero esto es sólo porque la densidad normal tiene conjuntos de nivel que son elipses cuya forma depende directamente de las varianzas y correlaciones de \({\bf X}\).
Si \({\bf X}=(X_{1},X_{2})'\) fuera un vector con densidad uniforme en el cuadrado \([-1/2,1/2]^2\), entonces las componentes \(X_{1}\) y \(X_{2}\) de \({\bf X}\) serían independientes entre sí y con distribución uniforme en el intervalo \([-1/2,1/2]\). En este caso es fácil comprobar que \({\boldsymbol\mu}= {\bf 0}\) y que la matriz de varianzas-covarianzas es \[ {\boldsymbol\Sigma} = 12 \left( \begin{array}{rr} 1 & 0 \\ 0 & 1 \end{array} \right). \] Por tanto, dado un punto \({\bf x}=(x_{1},x_{2})'\), \(d_M({\bf x},{\boldsymbol\mu})=12\|{\bf x}\|^2\). Es decir, la distancia de Mahalanobis en este caso es proporcional a la distancia euclídea.
¿Cómo es posible que el vector \(\mathbf X=(X_1,X_2)'\) tenga marginales normales y no sea normal bivariante?
Pongamos un ejemplo de distribución bivariante no gaussiana con marginales normales. Sea \(X_1\) una v.a. N(0,1) y sea \[ X_2 = \left\{ \begin{array}{rl} -X_1 & \mbox{si $-1\leq X_1\leq 1$;} \\ X_1 & \mbox{en otro caso.} \end{array} \right. \]
Calculando la función de distribución de \(X_2\), \(F_{X_{2}}(x_{2})=\mathbb P\{X_{2}\leq x_{2}\}\), es fácil ver que \(X_2\) sigue una distribución \(N(0,1)\): \[ {\mathbb P}\{X_{2}\leq x_{2}\} = \left\{ \begin{array}{ll} {\mathbb P} \{ X_{1} \leq x_{2} \} = \Phi(x_{2}) & \mbox{si $x_2< -1$} \\ {\mathbb P} \{ X_{1} < -1 \} + \mathbb P \{ x_2\leq X_{1}\leq 1\} = 1-\Phi(1)+\Phi(1)-\Phi(x_2) = \Phi(x_2) & \mbox{si $-1 \leq x_2\leq 1$} \\ {\mathbb P} 1-\{ X_{1} \geq x_{2} \}=1-(1-\Phi(x_2)) =\Phi(x_2) & \mbox{si $x_2> 1$} \end{array} \right. \] donde \(\Phi\) denota la función de distribución normal estándar.
Calculando la probabilidad de que \(X_1-X_2\) sea igual a 0 y viendo que no es 0, deducimos que la distribución conjunta de \(X_1\) y \(X_2\) no es normal bivariante, porque, si lo fuera, entonces \(X_1-X_2\) también debería ser normal univariante y en este caso \(\mathbb P\{X_1-X_2=0\}=0\): \[ \mathbb P\{X_1-X_2=0\} = \mathbb P \{X_1< -1\} + \mathbb P \{X_1> 1\} = 2(1-\Phi(1)) \neq 0. \]
Lo que no he calculado es \(\mbox{Cov}(X_1,X_2)\), pero creo que no será 0. Me faltaría pensar un ejemplo de dos marginales univariantes \(X_1\) y \(X_2\) que sean \(N(0,1)\) y con \(\mbox{Cov}(X_1,X_2)=0\), pero tal que \(\mathbf X=(X_1,X_2)'\) no sea normal bivariante. Si a alguien se le ocurre, que me lo diga.
En la descomposición espectral de una matriz de covarianzas \(\boldsymbol\Sigma = \bf C \bf D \bf C'\), siendo \(\bf C\) la matriz ortogonal cuyas columnas son los autovectores ortonormalizados de \(\boldsymbol\Sigma\) y \(\bf D\) la matriz diagonal con los correspondientes autovalores de \(\boldsymbol\Sigma\), ¿es \(\bf D\) la matriz de Jordan?.
Agradezco a Ernesto Girondo que me resolviera la pregunta. La respuesta es sí.
Queremos comprobar si la distribución de unas localizaciones en una ventana de observación es uniforme. Para ello dividimos la ventana superponiendo una malla que da lugar a una partición de rectángulos \(A_1,\ldots,A_k\), todos de la misma área. Definimos la v.a. cualitativa \(X\) = ‘’rectángulo en el que se observa la localización’‘, cuyo espacio muestral es \(A_1,\ldots,A_k\). También podemos definir la v.a. cuantitativa \(Y\) = “Número de localizaciones observadas en un rectángulo de la malla’’, cuyo espacio muestral es 0,1,2,… ¿Se usa un contraste de bondad de ajuste de la v.a. \(X\) a una uniforme o se contrasta si todas las \(Y\) sigue una distribución de Poisson?
Las dos opciones son correctas, pero depende de si tenemos los datos originales y, si no, de cómo están agrupados los datos. Si dividimos la ventana de observación de las localizaciones observadas mediante una retícula muy fina entonces seguramente funcione mejor el contraste de bondad de ajuste a una Poisson. Si la retícula tiene rectángulos muy amplios en cada uno de los cuales pueden caer muchas observaciones, entonces es mejor usar el contraste a una uniforme.
Quizá se entienda mejor con un ejemplo.
Genero \(n\) observaciones de una distribución uniforme en el cuadrado unidad.
n <- 100
Datos <- cbind(runif(n),runif(n))
par(pty="s")
plot(Datos,xlab="",ylab="",pch=20,xlim=c(0,1),ylim=c(0,1))
Separo los datos según una malla “grosera” de 3x3. La función
quadratcount
nos devuelve la frecuencia observada de cada
rectángulo.
library(spatstat)
Datosppp <- ppp(x=Datos[,1],y=Datos[,2],xrange=c(0,1),yrange=c(0,1))
Q <- quadratcount(Datosppp,nx=3,ny=3)
Q
## x
## y [0,0.333) [0.333,0.667) [0.667,1]
## [0.667,1] 9 8 13
## [0.333,0.667) 11 11 15
## [0,0.333) 10 10 13
plot(Q)
Contraste de bondad de ajuste de \(X\) a una uniforme: \[ H_{0}: \; {\mathbb P}\{X=A_j\} = \frac{1}{k}, \;\; j=1,\ldots,k. \]
vectorQ <- matrix(Q,nrow=1)
chisq.test(vectorQ)
##
## Chi-squared test for given probabilities
##
## data: vectorQ
## X-squared = 3.5, df = 8, p-value = 0.8992
Vemos que no podemos rechazar la hipótesis nula, que es la decisión acertada.
Contraste de bondad de ajuste de \(Y\) a una Poisson(\(\lambda\)): \[ H_{0}: \; Y \sim \mbox{Poisson}(\lambda) \;\; \mbox{para algún } \lambda>0. \] Estimamos el parámetro \(\lambda\) de la Poisson:
N <- length(vectorQ) # Número de rectángulos en la malla
lambdaEst = sum(vectorQ)/N # Número promedio de localizaciones por rectángulo de la malla
lambdaEst
## [1] 11.11111
Valores que ha tomado la variable \(Y\) en la muestra (es un valor de \(Y\) por cada rectángulo de la malla).
sort(vectorQ)
## [1] 8 9 10 10 11 11 13 13 15
Vemos que, aunque el espacio muestral de la Poisson es \({\mathbb N}\), la cantidad de valores observados de \(Y\) es pequeña, porque la partición es grosera. Sin embargo, nosotros tenemos que calcular las frecuencias observadas y esperadas (bajo \(H_0\)) de muchos valores del espacio muestral, tanto si los hemos observado como si no. Vamos a calcular las frecuencias de los valores 0,1,2,…,\(M\) = máximo valor observado de \(Y\).
M <- max(vectorQ)
# Frecuencias observadas y esperadas para los valores 0,1,2,...,M
obs <- rep(NA,M+1)
esp <- rep(NA,M+1)
for (i in 1:(M+1)){
obs[i] <- sum(vectorQ == (i-1))
esp[i] <- n * dpois(i-1,lambdaEst)
}
esp[M+1] <- n-sum(esp[1:M])
Estchi2 <- sum(((obs-esp)^2)/esp)
Estchi2
## [1] 83.3726
pvalor <- pchisq(Estchi2,df=(M+1)-1-1,lower.tail=F)
pvalor
## [1] 6.669131e-12
Vemos que el p-valor es muy pequeño, por lo que rechazamos la hipótesis nula, que es una decisión incorrecta.
Ahora repito el proceso anterior con una malla “fina” de 10x10. El contraste de bondad de ajuste de \(X\) a la uniforme ni siquiera es válido, porque las frecuencias observadas son muchas menores que 5. Por el contrario, el contraste de bondad de ajuste de \(Y\) a la Poisson funciona estupendamente (en general): el p-valor no es significativo.
Q <- quadratcount(Datosppp,nx=10,ny=10)
Q
## x
## y [0,0.1) [0.1,0.2) [0.2,0.3) [0.3,0.4) [0.4,0.5) [0.5,0.6) [0.6,0.7)
## [0.9,1] 1 0 1 3 1 1 2
## [0.8,0.9) 0 1 0 0 1 0 0
## [0.7,0.8) 2 1 1 1 0 0 1
## [0.6,0.7) 1 0 1 0 1 0 4
## [0.5,0.6) 1 2 2 0 3 0 1
## [0.4,0.5) 1 0 1 1 1 1 1
## [0.3,0.4) 0 0 1 1 1 0 1
## [0.2,0.3) 2 0 2 2 1 1 0
## [0.1,0.2) 0 0 0 1 0 0 1
## [0,0.1) 1 2 3 1 1 0 2
## x
## y [0.7,0.8) [0.8,0.9) [0.9,1]
## [0.9,1] 3 1 0
## [0.8,0.9) 1 4 0
## [0.7,0.8) 1 0 0
## [0.6,0.7) 0 0 2
## [0.5,0.6) 1 3 0
## [0.4,0.5) 2 1 3
## [0.3,0.4) 0 1 3
## [0.2,0.3) 4 1 0
## [0.1,0.2) 0 3 1
## [0,0.1) 0 0 2
vectorQ <- matrix(Q,nrow=1)
chisq.test(vectorQ)
## Warning in chisq.test(vectorQ): Chi-squared approximation may be incorrect
##
## Chi-squared test for given probabilities
##
## data: vectorQ
## X-squared = 108, df = 99, p-value = 0.252
N <- length(vectorQ) # Número de rectángulos en la malla
lambdaEst = sum(vectorQ)/N # Número promedio de localizaciones por rectángulo de la malla
lambdaEst
## [1] 1
sort(vectorQ)
## [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [75] 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4
M <- max(vectorQ)
# Frecuencias observadas y esperadas para los valores 0,1,2,...,M
obs <- rep(NA,M+1)
esp <- rep(NA,M+1)
for (i in 1:(M+1)){
obs[i] <- sum(vectorQ == (i-1))
esp[i] <- n * dpois(i-1,lambdaEst)
}
esp[M+1] <- n-sum(esp[1:M])
obs
## [1] 37 40 12 8 3
esp
## [1] 36.787944 36.787944 18.393972 6.131324 1.898816
Estchi2 <- sum(((obs-esp)^2)/esp)
Estchi2
## [1] 3.712438
pvalor <- pchisq(Estchi2,df=(M+1)-1-1,lower.tail=F)
pvalor
## [1] 0.2942367
La cuestión es que no entiendo qué muestra R
en
Deviance Residuals
. Por ejemplo, en el ejemplo de la
esclerosis múltiple nos muestra esto:
Deviance Residuals:
Min 1Q Median 3Q Max
-1.31250 -0.45006 -0.21017 0.01804 2.63976
En principio creía que mostraba información sobre la desviación en el individuo \(i\) (definida como \(D^2_i = -2[y_i \, \log \eta(x_i) + (1 - y_i) \, \log(1 - \eta(x_i))]\): su valor mínimo, su mediana, etc. Sin embargo, creo que estos valores son siempre positivos, por lo que no tiene sentido que su mediana sea un valor negativo.
Cada deviance residual del output de R
es la raíz
cuadrada de la desviación \(D^2_i\),
con un signo asociado: +1 si \(y_i=1\),
-1 si \(y_i=0\). Los valores de estos
residuos los puedes extraer con la función resid
:
reglog = glm(V6~V1+V2+V3+V4+V5,data=Datos,family="binomial")
summary(resid(reglog))
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.31250 -0.45006 -0.21017 -0.08950 0.01804 2.63976
El comando ClasFisher$scaling
¿da el vector
\({\mathbf w}_F\)? Por otro lado, ¿qué
significa el output que da summary(ClasFisher)
? Y por
último, no entiendo muy bien el procedimiento leave-one-out, ni para qué
sirve ni la forma en que lo implementas en R.
ClasFisher$scaling
no da el vector \({\mathbf w}_F\) exactamente, sino un vector
proporcional. Pero no sé decir por qué. Ni siquiera tiene necesariamente
longitud 1. De todas formas, para el procedimiento de clasificación no
importa la longitud por la expresión de la regla de decisión (que es una
desigualdad en la que en ambos lados se multiplica escalarmente por ese
vector de \({\mathbf w}_{F}\)).
¡Yo tampoco sé qué proporciona el output de
summary(ClasFisher)
! No lo hemos visto en clase y lo he
mirado ahora pero no sé qué da. En cualquier caso, no es nada
interesante. No merece la pena que perdamos el tiempo en averiguar en
qué consiste.
El procedimiento leave-one-out (LOO) es un caso particular de una
importante metodología práctica (computacional) que se
utiliza en Estadística para comprobar si un procedimiento (típicamente
de clasificación o regresión) funciona bien o no. LOO es un caso
particular de validación cruzada (o CV, por sus siglas en inglés,
cross-validation). CV consiste en un procedimiento iterativo que, en
cada iteración del bucle, separa la muestra en dos submuestras (una
grande y una pequeña). Con la grande construye un clasificador (si
nuestro problema estadístico es la clasificación) y lo utiliza para
clasificar la submuestra pequeña. Entonces, para la submuestra pequeña
se comparan las clases que el clasificador ha asignado a los individuos
con las clases a las que pertenecen realmente. Se guarda la información
del número de clasificaciones erróneas y el tamaño de la muestra pequeña
en cada iteración. Al terminar el procedimiento iterativo (un bucle for
en R
) se utiliza la información guardada para proporcionar
una estimación de la proporción o tasa de clasificaciones erróneas de
ese clasificador (en nuestro caso, el de Fisher o la regresión
logística).
Seguramente se entiende mejor en el caso del algoritmo LOO. En este caso, en la iteración \(i\) del bucle for se toma la observación \(i\)-ésima de la muestra. Esta observación es la submuestra pequeña y las restantes \(n-1\) observaciones son la submuestra grande. Con la submuestra grande construyo el clasificador y con éste clasifico la observación \(i\). Acumulo la clase en la \(i\)-ésima componente de un vector de tamaño \(n\) que he inicializado antes del bucle. Al finalizar el bucle tengo un vector de las poblaciones que el clasificador ha asignado a cada individuo de la muestra. Comparo la clase decidida por el clasificador con la clase a la que realmente pertenece cada individuo. Con esto determino la proporción de fallos (tasa de error) del clasificador. Si en lugar del usar LOO uso otro procedimiento de CV, esta tasa de error puede cambiar (esperemos que poco).
Respecto a cómo llevar a cabo LOO en la práctica para el clasificador
de Fisher con R
, conviene aprovechar el output de la
función lda
del paquete MASS
. El resultado
final es el mismo que en el algoritmo que expliqué arriba.
CV=TRUE
en la función
lda
.ClasFisher = lda(X,Y,CV=TRUE)
n = nrow(X) # Tamaño muestral
Y
) con la clase
asignada por el clasificador, que se obtiene con
ClasFisher$class
. Obtengo la proporción de clasificaciones
erróneas.sum(Y != ClasFisher$class)/n
¿Cómo divido la pantalla de gráficos de R en una cuadrícula?
Hay dos maneras. Supongamos que queremos dividirla en una cuadrícula
de 2 filas y tres columnas. Con layout
le tenemos que dar
una matriz con componentes de 1 a 6 en 2 filas y 3 columnas:
X = matrix(rnorm(20),ncol=2) # Datos
layout(matrix(c(1:6),2,3)) # Divido la pantalla de gráficos
plot(X,pch=21)
plot(X,pch=21)
plot(X,pch=21)
plot(X,pch=21)
plot(X,pch=21)
plot(X,pch=21)
layout(matrix(1)) # Deshago la división de la pantalla
Segunda manera: con mfrow
directamente le indico el
número de filas y columnas.
X = matrix(rnorm(20),ncol=2) # Datos
par(mfrow=c(2,3)) # Divido la pantalla de gráficos
plot(X,pch=21)
plot(X,pch=21)
plot(X,pch=21)
plot(X,pch=21)
plot(X,pch=21)
plot(X,pch=21)
par(mfrow=c(1,1)) # Deshago la división de la pantalla
Cuando escribo set.seed
luego siempre se genera
el mismo número aleatorio, aunque no vuelva a ejecutar
set.seed
.
No. Ha debido de ser que seguías ejecutando el set.seed
aunque no te dieras cuenta. Por ejemplo,
set.seed(20)
rnorm(1)
## [1] 1.162685
rnorm(1)
## [1] -0.5859245
rnorm(1)
## [1] 1.785465
Lo que sí es cierto es que la semilla fijada determina (secuencialmente) todos los números aleatorios que se generen a partir de ese momento. Es decir, si vuelvo a escribir las cuatro líneas de código de arriba volveré a obtener exactamente los mismo tres números aleatorios.
Normalmente al iniciar una nueva sesión de R
, el
programa internamente toma una semilla nueva basada, entre otras cosas,
en la hora. Para deshacer la orden de fijar la semilla se pueden usar
cualquiera de los dos comandos siguientes
rm(.Random.seed)
set.seed(NULL)
¿Cómo convierto una matriz de números en una tabla de datos?
Con la función data.frame
. Por ejemplo,
v = 1:6
M = matrix(v,nrow=3) # Matriz numérica 3x2
Datos = data.frame(M) # Convertimos la matriz M en una tabla de datos (objeto data.frame)
names(Datos) = c("a","b") # Le damos nombre a las columnas de la tabla