estadistica¶

2020-05-03
in estadistica
6 min read

ANOVA

library(tidyverse)
library(patchwork)

ANOVA refiere a "Analysis of Variance" en inglés y corresponde a una serie de procedimientos estadísiticos que permiten estudiar diferencias de medias poblacionales, basado en muestras.
Es una técnica muy difundida para comparar medias de 2 o más grupos. Específicamente queremos ver si todos los grupos comparten media o al menos uno difiere. En el caso más simple, de comparar dos medias, el resultado es equivalente al test t de comparación de medias por lo que ANOVA se considera una generalización de este.

El test de hipótesis sería: $\(H_0: \mu_1 = ... = \mu_k$\) $\(H_1: \text{las medias no son iguales}$\) ANOVA tiene también una serie de supuestos que hay que tener en cuenta.

Independencia de las observaciones.
Normalidad en los residuos. Podemos pensarlo como normalidad dentro de cada grupo, siendo el residuo la parte no explicada por la media del grupo. En muestras chicas puede ser problemático si no se cumple (reduce la potencia del test). Con muestras grandes debería cumplirse por Teorema Central del Límite.
Homocedasticidad. Se supone que cada grupo tiene misma varianza. Si la muestra no es muy chica ANOVA es bastante robusto con este supuesto, si no, hay alternativas no parámetricas por ejemplo.

Un ejemplo simulado

Generamos primero un set de datos donde la media de 3 grupos es distinta y vamos paso a paso con los cálculos.

set.seed(24)
df = data.frame(grupo1 = rnorm(n = 35, mean = 3, sd = 1),
                grupo2 = rnorm(n = 35, mean = 6, sd = 1),
                grupo3 = rnorm(n = 35, mean = 9, sd = 1))

Tenemos 3 grupos de 35 observaciones, cada uno proveniente de poblaciones con medias notoriamente distintas.
Veamos como resultaron las medias muestrales.

sapply(df, FUN = mean)

##   grupo1   grupo2   grupo3 
## 2.916989 5.910171 8.904245

Viédolo gráficamente vemos que son muy dispares y ANOVA debería captar estas diferencias.
En el segundo gráfico, la linea vertical represnta la media general del dataset

Veamos como resulta analizar esto con ANOVA.

Cálculos.

Obviamente existen paquetes estadísticos para realizar este análisis rapidamente pero iremos paso por paso.
La lógica es comparar la media de las poblaciones y para ello nos basamos en la varianza. Más precisamente en la descomposición de la varianza.

$\(SC_{Total} = SC_{Entre} + SC_{Dentro}$\) La suma de errores cuadrados de TODO el dataset se puede descomponer entre el desvío cuadrado de cada grupo frente a la media general (Entre) más el desvío cuadrado de cada observación respecto a su media grupal.

Siendo más intuitivos. Los suma de los desvíos cuadrados de cada observación respecto a la media general 5.91 pueden ser vistos como la diferencia entre medias grupales (qué tan lejos está cada pico del gráfico de la media grupal) más qué tan dispersos están los datos dentro de cada grupo.

Cuanto más grande sea la brecha entre la variabilidad entre grupos y la variabilidad al interior de los grupos, más probable es que las medias poblacionales sean distintas. Es decir, si la variabilidad total se explica más por la diferencia entre medias grupales que por la diferencia entre desviós al interior, entonces más evidencia en favor de distintas medias grupales. Si el ratio no es tan grande, entonces tenemos menos fuerza para afirmar tal cosa.

Para poder comparar correctamente, no se mira directamente $SC_{Entre}$ vs $SC_{Dentro}$ ya que estos dependen del tamaño de la muestra, sino que se los normaliza primero. $SC_{Entre}$ se normaliza por sus grados de libertad siendo k-1 (cantidad de grupos menos 1) y $SC_{Dentro}$ se normaliza con N-K (observaciones totales menos cantidad de grupos).

Detalle técnico

Como asumimos que los residuos son normales, entonces elevarlos al cuadrado nos devuelve una distribución Chi-Cuadrado. Las sumas de residuos al cuadrado son entones Chi-Cuadrado con los grados de libertad que mencionamos. Si dividimos dos distribuciones Chi-Cuadrado, normalizadas por sus grados de libertad, obtenemos una distribución F con grados de libertad equivalentes a los de ambas Chi-Cuadrado.

Ese estadístico F, que sigue la distribución recién mencionada sera nuestro estádistico para testear la Hipótesis.

\[F = \frac{\frac{SC_{Entre}}{K-1}}{\frac{SC_{Dentro}}{N-K}}\]

Donde: $\(SC_{Entre} = \sum_{i=1}^k{n_i (\bar{x}_i} - \bar{x})^2$\) $\(SC_{Dentro} = \sum_{i=1}^K\sum_{j=1}^{n_k}{(x_j - \bar{x}_i)^2}$\) Luego como en cualquier test de hipótesis, comparamos el estadístico F con la distribución teórica si la hipótesis nula fuera cierto y según el valor de alfa que hayamos elegido, rechazamos o no la hipótesis nula.
Para ilustrar, la dsitribución F tiene la siguiente forma con los grados de libertad de nuestro ejemplo.

Donde la región en rojo es el area de la curva posterior al 95% de la distribución. Si nuestro estadístico cae en la zona rojo podemos rechazar la hipótesis nula con alfa =0.05

Obtengamos los números con la función aov.

res = aov(valor ~ grupo, data = df_long)
summary(res)

##              Df Sum Sq Mean Sq F value Pr(>F)    
## grupo         2  627.3  313.66     326 <2e-16 ***
## Residuals   102   98.1    0.96                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La tabla que devuelve la función es justamente todo lo que fuimos viendo.
Sum Sq es la suma de desvíos cuadrados. La fila de grupo corresponde a Entre y Residuals corresponde a Dentro.
DF son los grados de libertad (K-1) y (N-K).
Mean Sq es la división de Sum Sq por sus grados de libertad. Serían el numerador y denominador del estadístico F.
F value es simplemente la división de los Mean Sq. Obtenemos un estadístico de 326(!). A partir de 3 aprox ya podíamos rechazar la hipótesis nula. El p-value (la última columna) es virtualmente 0.

Dado este resultado, podemos rechazar la hipótesis nula y asegurar con el 95% de confianza que las medias poblacionales no son iguales.
Era un caso medio extremo pero sirve de ejemplo.

Lo que no nos dice ANOVA es si todas son distintas o cuál es diferente al resto. Para ello hay que hacer estudios posteriores pero no entramos en detalles acá.

Un caso de medias iguales

res2 = aov(valor ~ grupo, data = df2_long)

summary(res2)

##              Df Sum Sq Mean Sq F value Pr(>F)
## grupo         2   0.00  0.0014   0.001  0.999
## Residuals   102  98.14  0.9621

Aquí es el otro extremo. El valor del estadístico F es casi 0, por lo tanto el p-value es casi 1. No hay evidencia para rechazar la hipótesis nula.

Conclusiones

ANOVA cómo método para comparar medias poblacionales es muy sencillo de aplicar y bastante robusto frente a inconsistencias en los supuestos. Permite dar una medida objetiva de si es posible o no rechazar la hipótesis nula, más allá de que uno pueda tener una primera impresión visual.

Para ir un paso más allá, ANOVA se puede relacionar directamnte con las regresiones lineales. Anova tal como lo presentamos es equivalente a correr una regresión donde la variable independiente es el grupo al que pertenece la observación. Las generalizaciones como ANCOVA, MANOVA, etc, también tienen su correlato en regresión. Esto sucede porque según el campo de estudio se eligieron caminos y convenciones distintos de análisis, llevando a distintas ramas que al final hacen lo mismo, pero genera confusión al intentar entender la estadística como un todo.

2020-04-25
in estadistica, R
5 min read

Distribucion Gamma

Origen y uso habitual

La distribución Gamma es continua y siempre positiva. Se parametriza con dos parámetros que deben ser positivos.

Lamentablemente no hay un consenso sobre cómo llamar a los parámetros y prevalecen dos formas, bastante similares pero con distinto origen.

La primera es con los parámetros shape $k$ y scale $\theta$.
La segunda es con los parámetros shape $\alpha$ y rate $\beta$.

La relación entre ambas es $k = \alpha$ y $\theta = 1/\beta$.
Obviamente cualquiera que se use va a resultar en los mismos resultados, pero hay que estar atento para lograr la parametrización adecuada.

Yo prefiero la segunda opción, donde el rate $\beta$ puede relacionarse al rate $\lambda$ de una poisson.

Por qué nos interesa relacionarla con la Poisson?

Uno de los usos habituales y que resulta fácil de entender es que si nos encontramos en un Proceso de Poisson (si no se recuerda ver ACA), suceden $\lambda$ eventos por período en promedio (la distribución poisson nos ayuda con eso), el tiempo (medido en períodos) entre un evento y otro se puede modelar con una exponencial con parámetro rate = $\lambda$ (el mismo de la poisson), y ahora adicionamos que el tiempo medido en períodos hasta que suceda el k-ésimo evento, se puede modelar con la distribución gamma, parametrizada con shape = $k$ (cantidad de eventos) y rate = $\lambda$ (el mismo de la poisson nuevamente.)

Por ejemplo:

Si podemos modelar la cantidad de veces que vamos por semana a comprar cerveza como una poisson con $\lambda$ = 2, es decir en promedio dos veces por semana, tendremos la siguiente distribución.

set.seed(1)
ggplot() +
  geom_bar(aes(x = rpois(10000,2)), fill = "lightgreen", color = "black") +
  xlab("") +
  ylab("") + 
  ggtitle("Simulación de 10000 variables Poisson con rate = 2") +
  scale_x_continuous(breaks = seq(0,9)) + 
  theme_minimal()

Donde la mayoría de las semanas vamos alrededor de 2 veces.

Podemos a su vez, modelar el tiempo entre cada evento con una exponencial. Se lo puede pensar como ¿cuánto tiempo falta para que vaya de nuevo a comprar cerveza a partir del momento en que estoy parado?

Usaremos una exponencial, pero para verlo en días, que sería lo apropiado, cambiamos $\lambda = 2$ por $\lambda = 2/7$, que sería el rate por día.

set.seed(2)
ggplot() +
  geom_histogram(aes(x = rexp(10000, 2/7)),bins = 100, fill ="darkgreen", color = "black") +
  xlab("") +
  ylab("") + 
  ggtitle("Simulación de 10000 variables Exponenciales con rate = 2/7") +
  scale_x_continuous(breaks = c(seq(0,10),15,20,25,30,35,40)) + 
  theme_minimal()

Vemos que en general faltan 1 o 2 días para que tengamos que ir de nuevo, aunque si tomamos el promedio veremos que es 3.5 dias, lo cuales lógico porque venimos de una poisson con rate de 2 veces por semana.

Si ahora queremos ver en cuantos días habremos ido 5 veces, podemos usar la distribución gamma, con $shape = 5$ (porque queremos ver el 5 evento), y $rate = 2/7$ porque es el rate diario.

set.seed(3)
ggplot() +
  geom_histogram(aes(x = rgamma(10000, shape = 5, rate =  2/7)),bins = 100, fill = "darkred", color ="black") +
  xlab("") +
  ylab("") + 
  ggtitle("Simulación de 10000 variables Gamma con shape = 5  y rate = 2/7") +
  scale_x_continuous(breaks = seq(0,60, by = 5)) + 
  theme_minimal()

Podemos ver que la distribución gamma es asimétrica en este caso y en general faltan unos 12-16 días. Si tomamos el promedio vemos que es de 17.5 días, lo cual tiene sentido ya que es 5 veces el tiempo promedio de espera, que era 3.5 días.

La media de una gamma puede calcularse rápidamente como $media = \frac{shape}{rate} = \frac{5}{2/7} = 17.5$

Hay una clara relación entre la distribución Exponencial y la Gamma.
Primero, como ya vimos, la exponencial modela el tiempo hasta el próximo evento en un proceso Poisson y la Gamma hasta el k-ésimo evento.
Podemos pensar a la distribución Gamma como la suma de K distribuciones Exponenciales con un mismo rate!

Yendo más allá, en realidad, la distribución Gamma es una familia de distribuciones, y la Exponencial no es más que un caso particular de la Gamma con k = 1.
Si Gamma nos permite saber el tiempo hasta el k-esimo evento, y la Exponencial es hasta el primer evento, entonces la Exponencial como caso particular de la Gamma parece obvio.

Un poco de formulas

La función de densidad de la distribución Gamma, utilizando $\alpha$ y $\beta$ es: $\(\frac{\beta^{\alpha}}{\Gamma (\alpha)}x^{\alpha - 1}e^{-\beta x}$\)

donde $\Gamma (\alpha)$ = $(\alpha - 1)!$

Cuidado con las distintas maneras de nombrar a los parámetros. En la literatura posiblemente vean $\alpha$ y $\beta$ cuando usen shape y rate. Aquí para el ejemplo de las cervezas reemplazamos $\alpha$ por $k$ y $\beta$ por $\lambda$ para relacionarlo con la distribución Poisson.

El equivalente con esa notación sería: $\(\frac{\lambda^{k}}{\Gamma (k)}x^{k - 1}e^{-\lambda x}$\)

Como dijimos si reemplazamos k por 1, obtenemos la función de densidad de la Exponencial.

Otras distribuciones que son casos particulares de la Gamma son la Erlang (es Gamma pero con valores discretos de K, la exponencial también es un caso particular de Erlang con k= 1) y la Chi-Cuadrado

Desestimar el significado de Shape y Scale

Ya sea que prefiramos shape y scale o shape y rate como parametrización, sus nombres en inglés llevan a pensar que la dsitribución varía su forma con shape y su escala con scale. Eso no es tan así y puede llevar más a confusiones que otra cosa. Ambos parámetros pueden afectar tanto la forma y escala, por eso me parece más sencillo pensarlo con $k$ Y $\lambda$, al menos en el ámbito de los procesos Poisson.
Igualmente, no olvidar que la distribcuión Gamma se puede usar para otros campos que no son situaciones de un Proceso de Poisson y esa interpretación de los parámetros puede ser poco apropiada.

k = c(2,9)
rate = c(1,2)
x= seq(0,10, by = 0.1)

df <- expand.grid(lambda = as.factor(rate), k = k, x = x) %>%
  mutate(gamma = dgamma(x,k, as.numeric(as.character(lambda))),
         combination = as.factor(paste0("k = ",k,", lambda = ",lambda)))

ggplot(data = df, aes(x = x, y =gamma, color = combination)) + 
  geom_line() + 
  #geom_point() + 
  ggtitle("Función de Densidad para distintos parámetros de una Distribución Gamma") +
  xlab("Períodos") + 
  ylab("") + 
  scale_color_brewer(palette = "Set1")

2020-04-18
in estadistica
5 min read

Proceso Poisson y distribucion exponencial

Basado en:
https://stats.stackexchange.com/questions/2092/relationship-between-poisson-and-exponential-distribution
https://towardsdatascience.com/the-poisson-distribution-and-poisson-process-explained-4e2cb17d459

Un proceso que sigue un proceso de Poisson es aquel que tiene un cantidad de eventos promedio cada determinada unidad de medida (generalmente tiempo). Estos sucesos deben ser independientes y aleatorios pero la cantidad debe estar centrada alrededor del promedio ya mencionado. Una última condición es que dos sucesos no pueden suceder en simultáneo.

Para ejemplificar, podemos suponer que la cantidad de pacientes que llegan a la guardia de un hospital sigue un proceso de Poisson.

Un proceso Poisson se puede "descomponer" en 2 conceptos.

Una disitribución Poisson (no proceso!) que provee la función de densidad de la cantidad de pacientes que ingresan al hospital. (Discreta)
Una distribución exponencial que modela el tiempo que transcurre entre cada paciente. (Continua)

Distribución Poisson

Esta distribución sirve para modelar la cantidad de pacientes. Digamos que ingresan en promedio $\lambda$ pacientes por hora.

La función de probabilidad (probability mass function en inglés ya que es discreta) es: $\(P(K) = e^{-\lambda}\frac{\lambda^k}{k!}$\)

Donde K es la cantidad de pacientes en una hora y $\lambda$ es como dijimos, la cantidad promedio que ingresa por hora. La función nos dice qué probabilidad hay de recibir K pacientes en una hora si en general llegan $\lambda$.

Distintos $\lambda$ devuelven obviamente distintas probabilidades para cada valor de K. Vemos que a medida que crece $\lambda$ la densidad se parece cada vez más a una normal.

lambda = c(1,2,4,6,10)
k = seq(0,max(lambda)*1.5)
df <- expand.grid(lambda = as.factor(lambda), k = k) %>%
  mutate(poisson = dpois(k, as.numeric(as.character(lambda))))

ggplot(data = df, aes(x = k, y =poisson, colour = lambda)) + 
  geom_line() + 
  geom_point() + 
  ggtitle("Probabilidad por intervalo (1 hora)") +
  xlab("Cantidad K de pacientes") + 
  ylab("")

$\lambda$ determina la forma de la dsitribución y como es de esperar, esta se centra alrededor del parámetro ya que es la cantidad promedio en el intervalo. Una propiedad interesante es que la varianza de la disitribución también es $\lambda$.
Los pacientes que llegan al hospital siguen una distribución $\sim P(\lambda)$.

Otra caracterísica es que uno puede escalar la distribución para cualquier intervalo. Es decir que si nuestro modelo era para cantidad de pacientes por hora, uno puede multiplicar $\lambda$ por 3 si quiere la distribución cada 3 horas, o dividir por 6 si la quiere cada 10 minutos por ejemplo. Esto es asi porque en realidad en la función de probabilidad de la distribución Poisson, $\lambda$ está multiplicado por el intervalo t,pero se simplifica y se lo asume 1. Luego uno puede derivar para el lapso que desee multiplicando.

Distribución Exponencial

Sabemos que llegan $\lambda$ pacientes por hora. Nos gustaría ahora saber la distribución para el tiempo de espera hasta que llegue el próximo paciente. Es decir, parados en el momento t, qué probabilidad hay de tener que esperar un minuto, 2, 10 hasta el próximo paciente?

El concepto es el siguiente. Si queremos saber la probabilidad de que haya que esperar al menos X minutos, es lo mismo que calcular la probabilidad de que la cantidad de pacientes en el momento t+x sea igual a la cantidad de momento t. Es decir, que no haya llegado nadie.
$N_{t+x} = N_{t}$.

Podemos calcular eso con la distribución Poisson. Es simplemente la probabilidad de 0 pacientes en un intervalo x.

$$P(N_{t+x} - N_t = 0) = e^{-\lambda x}\frac{\lambda x^0}{0!} = e^{-\lambda x} $$ La probabilidad de obtener un nuevo paciente en el intervalo x es 1 menos la probabilidad de no obtener ninguno, por lo tanto: $$ P(X_t \leq x) = 1 - P(N_{t+x} - N_t = 0) = 1 - e^{-\lambda x}$$ Este último resultado es la función de probabilidad acumulada de la exponencial. La probabilidad de tener que esperar X o menos minutos depende de $\lambda$ y x.

Con la derivada obtenemos la Función de distribución de la exponencial. $\(f(x, \lambda) = \lambda e^{-\lambda x}$\)

Si al hospital llegaran 10 personas por hora a la guardia, tendríamos la siguiente función de densidad de la exponencial. Tener en cuenta que $\lambda$ = 10 en la poisson (10 casos por unidad de tiempo (hora)), sin embargo, en la exponencial queremos medirlo en minutos, por lo que usamos un $\lambda$ de 10/60.

lambda = c(10)
t = seq(0,30, by = 1)
df2 <- expand.grid(lambda = as.factor(lambda), t = t) %>%
  mutate(exponencial = dexp(t, as.numeric(as.character(lambda))/60))

ggplot(data = df2, aes(x = t, y =exponencial)) + 
  geom_line(colour = "red") + 
  geom_point(colour = "red") + 
  ggtitle("Función de densidad de distribución exponencial: lambda = 10/60") +
  xlab("Minutos") +
  ylab("Densidad") +
  scale_x_continuous(breaks = seq(0,30, by = 5))

Por otro lado podemos ver la probabilidad acumulada hasta determinado minuto. Siguiendo el mismo ejemplo de 10 pacientes por hora (Poisson con $\lambda = 10$), la acumulada de la distribución exponencial tiene la siguiente forma.

Hay alrededor de 80% de chances que un paciente llegue en los próximos 10 minutos.

lambda = c(10)
t = seq(0,max(lambda)*3, by = 1)
df2 <- expand.grid(lambda = as.factor(lambda), t = t) %>%
  mutate(exponencial = pexp(t, as.numeric(as.character(lambda))/60))

ggplot(data = df2, aes(x = t, y =exponencial)) + 
  geom_line(colour = "blue") + 
  geom_point() + 
  ggtitle("Probabilidad de tener que esperar X minutos o menos") +
  xlab("Minutos") + 
  ylab("") +
  scale_x_continuous(breaks = seq(0,30, by = 5))

Por último, si lo quieren ver al revés. Podemos ver la probabilidad de tener que esperar al menos X minutos para que llegue el próximo.

Hay alrededor de 5% de chances de tener que esperar 20 minutos hasta el próximo paciente

lambda = c(10)
t = seq(0,max(lambda)*3, by = 1)
df2 <- expand.grid(lambda = as.factor(lambda), t = t) %>%
  mutate(exponencial = pexp(t, as.numeric(as.character(lambda))/60, lower.tail = FALSE))

ggplot(data = df2, aes(x = t, y =exponencial)) + 
  geom_line(colour = "darkgreen") + 
  geom_point() + 
  ggtitle("Probabilidad de tener que esperar al menos X minutos") +
  xlab("Minutos") + 
  ylab("") +
  scale_x_continuous(breaks = seq(0,30, by = 5)) + 
  scale_y_continuous(breaks = seq(0,1, by = 0.1))

2019-11-11
in algebra, estadistica, matematica
7 min read

Distintas Distancias

library(tidyverse)

Si tenemos un espacio euclideo, es decir, una linea, un plano o un hiperplano, que son los espacios típicos de la geometría clásica, podemos calcular la distancia entre dos puntos que se hayen en él.

Es decir, cuál es la distancia entre los puntos A (1,1) y B (1,0) en un plano? Empecemos pensando en los casos donde todos los valores del vector son numéricos.

A = c(0,0,1,1)
B = c(0,0,1,0)
recta = c(1,1,1,0)
df = as.data.frame(matrix(data = c(A,B, recta),
                          nrow = 4, 
                          ncol = 4,
                          byrow = TRUE )) %>%
  rename( x0 = V1,
          y0 = V2,
          x1 = V3,
          y1 = V4)


ggplot(data=df[1:2,], aes(x=x0, y=y0)) + 
  geom_segment(aes(xend=x1, yend=y1),
               arrow = arrow(length = unit(0.3,"cm"))) + 
  geom_point( aes(x = x1, y = y1), 
              color = "red", size = 2)

Distancia Euclideana

La métrica más habitual que se utiliza es la distancia euclideana, que consiste en la recta que une ambos puntos.

ggplot(data=df[1:2,], aes(x=x0, y=y0))+
  geom_segment(aes(xend=x1, yend=y1),
               arrow = arrow(length = unit(0.3,"cm"))) + 
  geom_point(aes(x = x1, y = y1), 
               color = "red", size = 2) + 
  geom_segment(data = df[3,], 
               aes(xend=x1, yend=y1),
               color = "blue", 
               arrow = arrow(length = unit(0.3,"cm")))

Esta distancia se calcula con:
$\(d(A,B) = d(B,A) = \sqrt{(A_1 - B_1)^2 + (A_2 - B_2)^2 + ... + (A_n - B_n)^2} = \sqrt{\sum_{i=1}^n (A_i - B_i)^2}$\)

Como se ve en la imagen, los puntos A y B pueden verse como vectores que inician en el origen (0,0). La distancia euclidea es a su vez la distancia entre sus puntas, que a su vez puede pensarse como un vector de desplazamiento (de A a B por ejemplo).

En este caso la distancia euclidea es: $\(d(A,B) = \sqrt{ (1-1)^2 + (1 - 0)^2} = 1$\) Y que es algo visible en el gráfico.

De manera más general, podemos definir toda una familia de distancias en el espacio euclideo. Las distancias de Minkowsky.

La distancia Minkowsky de orden p es: $\(d(A,B) = d(B,A) = \Bigg({\sum_{i=1}^n |A_i - B_i|^p}\Bigg)^{1/p}$\)

Vemos que si p = 2, entonces la distancia de Minkowsky no es otra que la distancia euclideana.

Distancia de Manhattan

Otro valor que suele tomarse para p es p = 1, y eso corresponde a la distancia de Manhattan.

Esta distancia se calcula con:
$\(d(A,B) = d(B,A) = |A_1 - B_1| + |A_2 - B_2| + ... + |A_n - B_n| =\sum_{i=1}^n |A_i - B_i|$\)

Es básicamente la suma de las diferencias absolutas entre las distintas dimensiones de los vectores.

Luce asi.

A = c(0,0,3,3)
B = c(0,0,2,1)
recta = c(2,1,3,3)
manhattan1 = c(2,1,3,1)
manhattan2 = c(3,1,3,3)

df = as.data.frame(matrix(data = c(A,B, recta, manhattan1, manhattan2),
                          nrow = 5, 
                          ncol = 4,
                          byrow = TRUE )) %>%
  rename( x0 = V1,
          y0 = V2,
          x1 = V3,
          y1 = V4)


ggplot(data=df[1:2,], aes(x=x0, y=y0)) + 
    geom_point( aes(x = x1, y = y1), 
              color = "red", size = 2) + 
      geom_segment(data = df[3,], 
               aes(xend=x1, yend=y1, color = "blue"),
               #color = "blue", 
               arrow = arrow(length = unit(0.3,"cm"))) +
      geom_segment(data = df[4,], 
               aes(xend=x1, yend=y1,  color = "green",),
               #color = "green", 
               arrow = arrow(length = unit(0.3,"cm"))) + 
      geom_segment(data = df[5,], 
               aes(xend=x1, yend=y1),
               color = "green", 
               arrow = arrow(length = unit(0.3,"cm"))) +
      scale_colour_manual(name = 'the colour', 
         values =c('blue'='blue','green'='green'),
         labels = c('Euclideana','Manhattan'))

Vemos como el valor abosluto imposibilita ir en dirección diagonal. Lo que se logra es medir la distancia como si hubiera una grilla como la del gráfico. Su nombre proviene de su utilización para medir distancias al interior de una ciudad (uno no puede cruzar las manzanas por el medio!).

Para saber cual conviene utilizar hay que pensar en el problema en cuestión.

Ya sea medir distancias en ciudades o donde haya restricciones de ese tipo puede que Manhattan sea más apropiado.
Por otra parte al no elevar al cuadrado le da menos pesos a las grandes distancias o mismo outliers por lo que puede ser otro motivo válido.
Por último, algunos trabajos argumentan que es más adecuada en problema de alta dimensionalidad (o mismo valores menores a 1 en el exponente de la formula de Minkowsky)

Similaridad coseno

La similaridad coseno se utiliza cuando se quiere ver la similitud "angular" entre dos observaciones y no la distancia en el plano. Es decir, vemos la dirección pero no la magnitud

A = c(0,0,1,1)
B = c(0,0,2,2)
C = c(0,0,5,0)

df = as.data.frame(matrix(data = c(A,B,C),
                          nrow = 3, 
                          ncol = 4,
                          byrow = TRUE )) %>%
  rename( x0 = V1,
          y0 = V2,
          x1 = V3,
          y1 = V4)


ggplot(data=df[1:3,], aes(x=x0, y=y0 )) + 
  geom_segment(aes(xend=x1, yend=y1),
               arrow = arrow(length = unit(0.3,"cm"))) + 
  geom_point( aes(x = x1, y = y1), 
              color = "red", size = 2) + 
  geom_text(aes(x=x1, y = y1, label = c("A","B","C")),
            vjust = -0.5)

Si hicieramos la distancia euclideando entre A y B obtendriamos el valor de la distancia en el plano, sin embargo podemos ver que se encuentran sobre la misma recta y por lo tanto su dirección es la misma. La similaridad coseno mide el ángulo entre dos puntos. En este caso el ángulo entre A y B es 0, y por ende su similaridad coseno es 1. Ambas tendrían la misma similaridad con cualquier otro punto de la misma recta, por más alejado que esté. Respecto a C, tanto A y B tiene comparten el ángulo por lo tanto la similaridad coseno entre A y C será la misma que entre B y C.

cosA = c(1,1)
cosB = c(2,2)
cosC = c(5,0)

# Similaridad coseno entre A y B
lsa::cosine(cosA, cosB)[[1]]

## [1] 1

# Similaridad coseno entre A y C
lsa::cosine(cosA, cosC)[[1]]

## [1] 0.7071068

# Similaridad coseno entre B y C
lsa::cosine(cosC, cosB)[[1]]

## [1] 0.7071068

Hay que tener en cuenta el contexto de nuestro problema para decidir qué medida de distancia usar. Por ejemplo la similaridad coseno se usa de manera estándar en análisis de texto (text mining).

Distancia de Mahalanobis

La distancia de Mahalanobis tiene la particularidad que mide la distancia entre un punto (P) y una distribución de datos (D). Si tenemos una nube de puntos correspondiente a una distribución D, cuanto más cerca esté P del centro de masa (o "promedio") más cerca se encuetran P y D. Intuitivamente sirve para pensar si P puede pertenecer a D o no.
Dado que la nube de puntos no tiene por qué ser una esfera (donde cada dirección tiene la misma cantidad de puntos), hay que tener en cuenta cómo se dispersan los puntos alrededor del centro de masa.

No es lo mismo,

esfera = as.data.frame(MASS::mvrnorm(1000, mu = c(3,3), 
                                     Sigma = matrix(c(1,0,0,1),
                                                    nrow = 2,
                                                    ncol = 2)))

ggplot(data = esfera, aes(x = V1, y = V2)) + 
  geom_point() + 
  geom_point(data = as.data.frame(matrix(c(6,2),ncol = 2)), 
             aes(x = V1, y = V2), color = "red") + 
  geom_text(data = as.data.frame(matrix(c(6,2),ncol = 2)),
            aes(x = V1, y = V2,label = "P"),
            vjust = 1.5, color = "blue") +
  labs(title = "Distribución esférica")

que,

elipse = as.data.frame(MASS::mvrnorm(1000, mu = c(3,3), 
                                     Sigma = matrix(c(1,0.6,0.6,1),
                                                    nrow = 2,
                                                    ncol = 2)))



ggplot(data = elipse, aes(x = V1, y = V2)) + 
  geom_point() + 
  geom_point(data = as.data.frame(matrix(c(6,2),ncol = 2)), 
             aes(x = V1, y = V2), color = "red") + 
  geom_text(data = as.data.frame(matrix(c(6,2),ncol = 2)),
            aes(x = V1, y = V2,label = "P"),
            vjust = 1.5, color = "blue") +
  labs(title = "Distribución Elíptica")

Los centros de masa son los mismos y lo único que cambia es la matriz de variancias y covarianzas (o como se correlacionan las variables). La distancia de P al centro es la misma, pero está claro que en el caso esférico P se encuentra más cerca de la distribución que en el caso elíptico.

Mahalanobis tiene en cuenta este aspecto ya que involucra la matriz de varianzas y covarianzas.

La distancia entre el punto x y la distribución con vector de medias $\vec{\mu}$ y matriz de covarianzas S es: $$ D_M(\vec{x}) = \sqrt{(\vec{x} - \vec{\mu})^TS)$$}(\vec{x} - \vec{\mu})

Tanto el vector $\vec{x}$ como la distribución pueden ser multivariadas (como se ve en los gráficos de arriba).

Tener en cuenta que si tenemos dos puntos provenientes de la misma distribución, podemos usar la distancia de Mahalanobis como una medida de disimilaridad: $$ D_M(\vec{x},\vec{y}) = \sqrt{(\vec{x} - \vec{y})^TS)$$ Veamos por ejemplo como queda la distancia de P a las distribuciones esféricas y elípticas graficadas.}(\vec{x} - \vec{y})

# Caso Esférico

mahalanobis(x = c(6,2), 
            center = c(3,3), 
            cov = matrix(c(1,0,0,1),
                         nrow = 2,
                         ncol = 2))

## [1] 10

# Caso Elíptico
mahalanobis(x = c(6,2), 
            center = c(3,3), 
            cov = matrix(c(1,0.6,0.6,1),
                         nrow = 2,
                         ncol = 2))

## [1] 21.25

Queda claro que P es más cercano a la distribución esférica que a la elíptica.

2019-10-31
in estadistica, R
7 min read

Maxima Verosimilitud y estimacion bayesiana

Distribucion prior

A falta de una buena traducción usamos este término.

Supongamos que se toman muestras aleatorias de una distribucion con pdf (funcion de densidad de probabilidad) $f(x|\theta)$. Por ejemplo podrían provenir de una normal con media = $\mu$ y varianza = 4.
Nosotros no sabemos el valor de $\mu$ pero podemos tener una idea de qué valores puede tomar y tener en mente una distribución prior de este parámetro $\epsilon(\theta)$. Para el ejemplo sería $\epsilon(\mu)$. Podemos suponer que $\mu$ se distribuye como una uniforme (0,1) por decir algo.
El concepto radica en tener una distribución prior para los parámetros de la distribución de la cual tomamos muestras aleatorias.

Distribución Posterior

Volviendo a nuestra muestra $X_1...X_n$ proveniente de $f(x|\theta)$, podemos decir, debido a que son observaciones aleatorias e independientes que su distribución conjunta es $f_n(x_1...X_n|\theta) = f(x_1|\theta)...f(x_n|\theta)$, que lo podemos escribir como $f_n(x|\theta)$.
Dado que suponemos que $\theta$ proviene de una distribución $\epsilon(\theta)$, la pdf conjunta $f_n(x|\theta)$ tiene que ser vista como la pdf conjunta condicional de$X_1...X_n$ para un valor particular de $\theta$.
Multiplicando la pdf conjunta condicional por la pdf $\epsilon(\theta)$ obtenemos la (n+1) distribución conjunta de $X_1...X_n$ y $\theta$ bajo la forma $f_n(x|\theta)\epsilon(\theta)$. Sería la pdf de encontrar en simultáneo determinados valores para x y $\theta$. La probabilidad conjunta marginal de $X_1...X_n$ se encuentra integrando la pdf conjunta con $\theta$ para todos los valores de $\theta$. Sería la probabilidad marginal de encontrar determinados valores de x (sabiendo la distribución de $\theta$ pero sin saber el valor puntual que toma).

$g_n(x) = \int_\Omega f_n(x|\theta)\epsilon(\theta) d\theta$

Por teorema de Bayes tenemos que la distribución posterior de $\theta$, es decir, dados los x es: $\(\epsilon(\theta|x) = \frac{f_n(x|\theta)\epsilon(\theta)}{g_n(x)} \text{ para } \theta \in \Omega$\) Se dice que la distribución prior $\epsilon(\theta)$ representa la verosimilitud, antes de ver los valores de $X_1...X_n$, de que el verdadero valor de $\theta$ se encuentre en cada una de las regiones de $\Omega$ y que la pdf de la distribución posterior $\epsilon(\theta|x)$ representa la verosimilitud después que los valores $X_1 = x_1,...,X_n = x_n$ hayan sido observados.

## La funcion de Versoimilitud

El denominador de la distribución posterior es básicamente la integral del numerador para todos los posibles valores de $\theta$. Depende de los valores observados $X_1...X_n$ pero no de $\theta$, por lo que puede considerarse constante en este contexto.
Dado que es una constante podemos quitarla de la distribución posterior que vimos y decir que $\(\epsilon(\theta|x) \propto f_n(x|\theta)\epsilon(\theta)$\)

Cuando se ve $f_n(x|\theta)$ para una muestra aleatoria como función de $\theta$, se la suele llamar función de verosimilitud. En inglés: Likelihood function.

Juntando estos términos podemos decir que la pdf posterior de $\theta$ es proporcional al producto de la función de verosimilitud y la pdf prior de $\theta$.

La idea de ver esta relación de proporcionalidad es para poder calcular la pdf posterior evitando calcular la integral del denomiador $g_n(x)$. Si el numerador tiene la forma de alguna de las distribuciones conocidad (normal, beta, gamma, uniforme, etc) es posible calcular fácilmente el factor constante por el cual multiplicar esa pdf para llegar a la posterior.

Distribuciones prior Conjugadas

Este concepto refiere a que ciertas distribuciones son particularmente útiles para los cálculos cuando las variables aleatorias observadas provienen de alguna distribución específica.
Es decir que según la distribución de la que provienen las X puede que haya alguna distribución conjugada tal que al asumirla para la pdf prior $\epsilon(\theta)$ ya sabemos que la distribución posterior también será de esa familia.

Un ejemplo ilustrador:
Supongamos que tomamos observaciones $X_1...X_n$ de una distribución Bernoulli de la cual no sabemos el parámetro $\theta$ (que debe estar entre 0 y 1). Supongamos además que la pdf prior de $\theta$ es una distribución beta con algúnos parámetros dados $\alpha \text{ y } \beta$. En este caso sabemos que por ser un caso de distribución conjugada, la pdf posterior de $\theta$ dado $X = x_i (i = 1,...,n)$ es a su vez una distribución beta con parámetros $\alpha + \sum_{i=1}^n x_i \text{ y } \beta + n - \sum_{i=1}^n x_i$.

Según la distribución de la que provengan las observaciones hay distintas distribuciones conjugadas que son las más convenientes para ese caso.

Estimación de parámetros

La idea es estimar algún parámetro de la distribución de la cual se obtienen los datos observados. El valor estimado del parámetro va a depender de dos cosas:

Del estimador que hayamos elegido (es decir, la función de los datos elegida)
De la muestra. El valor estimado va a depender de los datos aleatorios que tengamos de la distribución.

Como el estimador depende de la muestra podemos verlo a su vez como una variable aleatoria.

Función de pérdida

Lo que queremos de un estimador es que devuelva un valor estimado "a" para el parámetro lo más cercano posible al verdadero valor de $\theta$. La función de pérdida es una función que cuantifica esto. $$ L(\theta,a)$$ Hay algunas funciones habituales pero pueden adecuarse según el problema.
Podemos decir que en general lo que se busca es encontrar una estimación para la cual la esperanza de la pérdida sea un mínimo.

Estimador bayesiano

Si tenemos una muestra aleatoria y una pdf posterior para $\theta$ entonces el valor esperado de la pérdida para cualquier valor estimado "a" es: $\(E[L(\theta,a)|x] = \int_\Omega L(\theta,a)\epsilon(\theta,x)d\theta$\)

Lo que buscamos es encontrar un valor de a cuya pérdida esperada sea mínima. La función que genera un valor de a mínimo para cada posible valor de X será un estimador de $\theta$ y en particular se llamará Estimador Bayesiano.
El estimador bayesiano, que minimiza la pérdida esperada para cualquier set de datos X, va a depender de la función de pérdida que elijamos y de la pdf prior que se elija para $\theta$.

Por ejemplo,para la función de pérdida más utilizada, que es la de error cuadrático $\(L(\theta,a) = (\theta -a)^2$\) está demostrado que la pérdida es mínima cuando $a$ es la media de la distribución posterior $E(\theta|x)$.

Dijimos que el valor del estimador bayesiano va a depender de la distribución prior elegida. Esto es cierto, pero hay que tener en cuenta que para muestras grandes las diferencias empiezan a achicarse y los estimadores bayesianos provenientes de distintos priors empiezan a converger en la mayoría de los casos.

Estimadores de Máxima Verosimilitud

Los estimadores de máxima verosimilitud (MLE) son muy comunmente usados para estimar parámetros desconocidos ya que más allá de la discusión casi filosófica de "bayesianos vs frecuentistas", sirven para estimar sin tener que definir una función de pérdida ni una distribución prior para los parámetros. Esto último es importante ya que para casos donde se necesita estimar un vector de parámetros, la distribución prior debe ser una multivariada que englobe a todos y eleva la complejidad del proceso bayesiano ampliamente.
Para muestras chicas MLE suele hacer un trabajo decente y para muestras grandes suele ser excelente por lo que se llega a resultados muy similares a través de un proceso más directo y más sencillo.

Para estimar mediante MLE lo único que necesitamos es la función de verosimilitud ya definida. $\(f_n(x_1...X_n|\theta)$\) Luego lo único que se hace es buscar el parámetro $\hat \theta$ (estimado) que maximice esa función. Básicamente es buscar qué parámetro hace que la probabilidad conjunta de obtener esos valores de X sea máxima? Ese es nuestro MLE.

Para la gran mayoría de los casos esta metodología funciona pero hay que tener en cuenta que es posible que para algunos problemas no haya un máximo para la función de verosimilitud o que haya más de un punto, en cuyo caso hay que elegir alguno de ellos.

MLE en Bernoulli

Supongamos que tomamos observaciones $X_1...X_n$ de una distribución Bernoulli de la cual no sabemos el parámetro $\theta$ (que debe estar entre 0 y 1).

Para cualquier vector de observaciones $X_1...X_n$ la función de verosimilitud es: $$ f_n(x|\theta) = \prod_{i = 1}^n \theta^{{x_i}(1-\theta)}$$ El valor de $\theta$ que maximice la función de verosimilitud es el mismo valor que maximiza $log f_n(x|\theta)$, por lo que es conveniente encontrar tal valor buscando que maximice: $\(L(\theta) = log f_n(x|\theta) = \sum_{i=1}^n[x_i log \theta + (1 - x_i) log(1-\theta)] = (\sum_{i=1}^nx_i)log \theta + (n-\sum_{i=1}^n x_i) log (1-\theta)$\)

Si derivamos $dL(\theta) / d\theta$ e igualamos a 0, resolviendo esa ecuando para $\theta$ encontramos que $\theta = \bar x_n$.
Este valor maximiza el logaritmo de la función de verosimilitud y por ende también de la función de verosimilitud en sí misma. Por lo tanto el MLE de $\theta$ es $\hat \theta = \bar X_n$

# Generamos 100 observaciones de una Bernoulli
set.seed(150)
data = rbinom(100, 1, prob = 0.723)

# Calculamos su promedio, que ya sabemos es la mejor estimación para p dados los datos
mean(data)

## [1] 0.68

# Definimos función de verosimilitud
# Es la pdf de una Bernoulli para cada observación y sumamos sus logaritmos (en negativo porque 
# el optimizador minimiza en vez de maximizar)
LL = function( p){
   R = dbinom(x = data, size = 1, prob = p)

   -sum(log(R))  # Negativo porque log de probabilidades es <0.
 }

# Función que busca los parámetros que minimzan el negativo de la log verosimilitud
# Elegimos un valor inicial de p en el medio.
stats4::mle(LL, start = list(p = 0.5) )

## 
## Call:
## stats4::mle(minuslogl = LL, start = list(p = 0.5))
## 
## Coefficients:
##         p 
## 0.6799996

Vemos que la estimación por MLE es idéntica a la media. No corresponde con el verdadero valor del parámetro poblacional p debido a la muestra particular que fue seleccionada.

Algunos comentarios finales:

En algunos casos no es posible encontrar la solución óptima si no es por métodos numéricos.
Cuando $n \to \infty$ MLE converge en probabilidad al verdadero $\theta$. Por ende cuando $n \to \infty$ el estimador bayesiano (que cumple la misma propiedad) y MLE serán muy parecidos entre sí y al verdadero $\theta$.
MLE solo depende de las observaciones y no de cómo y en qué orden fueron recolectadas.

2019-10-13
in estadistica, R
3 min read

Teorema Central del Limite

El teorema central del límite (TCL) es fundamental en el desarrollo de la estadística y ha obtenido distintas variantes a lo largo de la historia. Veremos dos de las versiones más conocidas.

Teorema Central del Límite para Media Muestral (Lindeberg - Lévy)

Si las varaibles $X_1 ... X_n$ forman una muestra aleatoria de tamaño n proveniente de una distribución con media $\mu$ y varianza $\sigma^2$ (0 < $\sigma^2$ < $\infty$), entonces para cualquier número fijo x. $$ \lim_{n\to \infty} Pr\Big[\frac{n^{½}(\bar X_n - \mu)}{\sigma} <= x\Big] = \Phi (x)$$

Donde $\Phi (x)$ es la función de distribución de una Normal Estándar.

El por qué de la convergencia del teorema no será probado acá pero no es díficil de encontrar. Por ejemplo ACÁ

Básicamente lo que dice el teorema, es que tomando una muestra grande de una población con media $\mu$ y varianza $\sigma^2$ definidas, entonces $\frac{n^{1/2}(\bar X_n - \mu)}{\sigma}$ va a tender a una normal estándar. Como consecuencia de eso podemos decir que $\bar X_n$ va a distribuirse aproximandamete como $N(\mu, \sigma^2/n)$.

El TCL nos dice cómo se distribuye la media muestral si tenemos una muestra grande.

Análogamente, también podemos decir que $\sum_{i=1}^n X_i$ va a ser aproximadamente una normal $N(n\mu, n\sigma^2)$

Ejemplo. Lanzar una moneda

Si lanzamos una moneda 900 veces. Cuál es la probabilidad de obtener más de 495 caras?

$X_i$ = 1 si sale cara en el lanzamiento i, y 0 si sale cruz.
E($X_i$) = ½ y Var($X_i$) = ¼. Esto se deduce de ser un experimento con distribución Bernouilli.

Para llevarlo a los términos del TCL, tenemos una muestra de tamaño n = 900, con $\mu$ = ½ y $\sigma^2$ = ¼.

Por TCL tenemos que la distribución de la suma del número total de caras $\sum_{i=1}^{900} X_i$ se distribuye aproximádamente como una normal con media = 900 * (½) = 450, varianza = 900 * (¼) = 225 y desvío estándar 225^(½) = 15.

Por lo tanto la variable $Z = \frac{H - 450}{15}$ se dsitribuye aproximadamente como una normal estándar. $\(Pr( H > 495) = Pr(\frac{H - 450}{15} > \frac{495 - 450}{15}) = Pr(Z>3) = 1 - \Phi(3) = 0.0013$\)

Podemos comparar contra el resultado que obtenemos al hacer el mismo ejercicio pero mirando directamente la distribución binomial (que es la que realmente genera el proceso de datos)

pbinom(495,900, 0.5, lower.tail = FALSE)

## [1] 0.001200108

Vemos que los resultados son muy similares.

Teorema Central del Límite para Suma de Variables Aleatorias Independientes (Liapunov)

Este TCL aplica a una secuencia de variables aleatorias independientes pero que no necesariamente tienen que provenir de una misma distribución. Todas deben tener una media y varianza definidas.

La variable $\(Y_n = \frac{\sum_{i=1}^n X_i - \sum_{i=1}^2 \mu_i}{(\sum_{i=1}^n\sigma_i^2)^{1/2}}$\)

Entonces $E(Y_n) = 0$ y $Var(Y_n)$ = 1

Siendo un poco más precisos veamos el teorema:

Suponiendo que las variables aleatorias $X_1. X_2, ...$ son independientes y que $E(|X_i - \mu_i|^3) < \infty$ para 1,2,... Y suponidendo que $\(\lim_{n\to \infty} \frac{\sum_{i=1}^n E(|X_i - \mu_i|^3)}{(\sum_{i=1}^n \sigma^2_i)^{3/2}} = 0$\) Entonces, utilizando la variable Y definida previamente tenemos que $\(\lim_{n \to \infty} Pr(Y_n <= x) = \Phi(x)$\)

La interpretacaión del teorema es que si se cumple la condición de los 3eros momentos, entonces para valores grandes de n la distribución de $\sum_{i=1}^n X_i$ será aproximadamente normal con media $\sum_{i=1}^n \mu_i$ y varianza $\sum_{i=1}^n \sigma^2_i$.

Diferencias entre Lindeberg-Lévy y Liapunov

El teorema de Lindeberg-Lévy aplica para secuencias de variables aleatorias iid y solo requiere que la varianza de estas variables sea finita. En cambio el teorema de Liapunov aplica a secuencias de variables aleatorias independientes pero que no necesariamente provienen de una misma distribución. Requiere que el tercer momento de cada variable existe y cumple con la ecuación del teorema.

Efecto de TCL

Más allá de la utilidad para aproximar distribuciones y medias mediante una normal, el TCL aporta una posible explicación a por qué tantas variables se distribuyen aproximadamante como normales. Si muchas de las variables a medir pueden pensarse como sumas de otras variables es lógico que tiendan a verse como normales aunque las variables que se suman para darle origen provengan de distintas distribuciones.

2019-07-01
in estadistica, R
11 min read

Metodos de Resampleo - ISLR Capitulo 5

Los métodos de resampleo son indispensables en la estadística moderna ya que permiten ajustar modelos a diferentes muestras de un mismo set de entrenamiento con el fin de obtener mayor información del modelo. Por ejemplo puede ser de utilidad para ver la variabilidad del modelo en distintas muestras. Los dos métodos que se presentan en el capítulo son cross-validation y bootstrap. A grandes rasgos CV puede servir para estimar el test error de un modelo o para ajustar hiperparámetros del modelo como el nivel de flexbilidad. Por su parte bootstrap puede usarse para medir la precisión de un parámetro estimado mediante un modelo estadístico.

Cross-Validation

De los modelos que uno entrena es de sumo interés obtener el "test error" que sería el error promedio al predecir una nueva observación aplicando el modelo estadístico entrenado. Esto puede calcularse si tenemos un set de testeo puntualmente para ello pero no suele ser el caso lamentablemente. En general no se tienen tantos datos como para separar en sets como uno desearía y surgen distintas técnicas para estimar el test error basado solamente en los datos de entrenamiento. Algunas de estas técnicas estiman el test error ajustando el training error por algún factor mientras que otras separan el training set en subsets donde uno hace las veces de test set.

Validation Set

Un método muy utilizado es el del set de validación. Básicamente consiste en separar nuestro training set en dos sets, un "nuevo" training set y uno set de validación. Una práctica habitual es separar 70-30, pero va a depender de la cantidad de observaciones que tengan y no hay una regla estricta. Básicamente de los datos que tienen para entrenar el modelo separan una parte que va a ser el set de validación y entrenan el modelo con los datos restantes (70% por ejemplo). Luego se mide la precisión del modelo en el 30% restante (set de validación) que son datos que no fueron utilizados a la hora de ajustar el modelo. Si utilizamos el MSE (mean squared error) cómo medida del error, este va a ser nuestro test error estimado. Recordemos que es el MSE calculado con las predicciones en el set de validación. Por otra parte el set de validación también puede servir para ajustar algún hiperparámetro. Se pueden correr muchos modelos con distintos hiperparámetros y ver cuál tiene menor MSE en el set de validación.

Es un método muy sencillo y suele ser eficaz pero tiene dos potenciales problemas: * El MSE puede variar mucho según cómo dividieron las observaciones en training y validación. Otra segmentación puede dar resultados muy distintos. * No utilizás todos tus datos para ajustar el modelo y puede que eso lleve a sobreestimar el test error, que quizás sería menor si usaras todas las observaciones para entrenar el modelo.

Leave-One-Out Cross-Validation

LOOCV es un intento de solucionar los problemas del enfoque del set de validación [SPOILER: No es recomendado pero vale la pena conocerlo].
Este enfoque es llevar el set de validación al extremo. Lo que se hace es de nuevo separar nuestro training set en dos pero esta vez guardando una sola observación como validación y usando las n-1 restantes para entrenar el modelo. La idea es hacer esto n veces, dejando cada vez una observación distinta como validación. El test error estimado es el promedio de los MSE de cada predicción que se hizo de la observación de validación. Pensando en los problemas del set de validación, con LOOCV logramos usar casi todos los datos disponibles para entrenar el modelo (n-1 observaciones) por lo tanto deberíamos tener modelos con menos sesgo y no sobreestimar tanto el test error como con el enfoque de set de validación. Por otra parte con el set de validación podemos obtener resultados muy distintos según el azar de cómo dividamos nuestros datos. En LOOCV esto no pasa ya que todos nuestros modelos de entrenamiento van a ser practicamente iguales salvo por una observación cada vez. No hay azar en la división de training y validación. Enseguida vemos el mayor problema de este enfoque, que es computacional. Debemos ajustar n modelos y no solo uno. Dependiendo de nuestros datos y la complejidad de nuestro modelo esto puede demandar muchísimo tiempo/recursos.

K-Fold Cross-Validation

K-Fold CV es un punto intermedio entre ambos enfoques y es de lo más utilizado al día de hoy. Consiste en separar nuestros datos en K subsets de mismo tamaño. Se selecciona uno de esos K subsets y se lo deja como validación. Se entrena el modelo con los K-1 subsets y se predice en el de validación que separamos. Así K veces, dejando como validación cada vez uno subset distinto. El Test error estimado es el promedio de los MSE en cada caso. Se puede ver fácilmente que si K = n, entonces estaríamos en LOOCV. Los valores típicos de K suelen ser 5 o 10, y por ende es muchísimo menos costoso que LOOCV. Al separar en "solo" 10 subsets cada set de validación puede tener cierta variabilidad en el MSE respecto a otros pero esta va a ser menor que en el enfoque de set de validación. En el libro se muestran unos gráficos para data simulada donde se ve que LOOCV y K-Fold tienen comportamiento muy similar y según el caso pueden sobreestimar o subestimar el verdadero test error (depende el problema y la flexibilidad elegida). Como mencionamos para el set de validación, K-fold también puede ser utilizado para ajustar algún hiperparámetro del modelo como el nivel de fleixibilidad. En este caso lo que nos interesa es encontrar el valor mínimo del MSE entre los distintos posibles valores del hiperparámetro para decidir cual es el mejor posible pero el valor puntual del MSE o su precisión no nos interesa tanto.

Trade-Off entre sesgo y varianza en K-Fold Cross-Validation

Otro punto muy importante de K-Fold, además de que requiere menos intensidad computacional que LOOCV, es que suele dar estimaciones más precisas del test error que LOOCV, y esto tiene que ver por el tradeoff entre sesgo y varianza.

Vimos antes que LOOCV debería ser el estimar más insesgado del test error ya que utiliza casi todas las observaciones de entrenamiento cada vez sin embargo hay que ver que sucede con la varianza ya que es otro componente del MSE. (Más detalles en ISLR Cap 2).
Resulta que LOOCV tiene mayor varianza que K-Fold CV siempre que K sea menor que n. Esto sucede porque en LOOCV lo que hacemos es promediar el resultado de n modelos cuyos datos de entrenamiento son casi idénticos (salvo por una observación) y por ende los resultados están en gran medida correlacionados positivamente.
Por otro lado al hacer K-Fold CV se promedian solo K resultados que están menos correlacionados entre sí ya que los datos de entrenamiento se solapan menos entre ellos. La clave acá es que el promedio de muchos valores altamente correlacionados tiene mayor varianza que el promedio de muchos valores que no están tan correlacionados. Dado este escenario se hicieron pruebas que llegaron a la conclusión empírica de que K=5 y k = 10 son valores que no suelen tener excesivo sesgo ni varianza.
Al parecer en los últimos años se empezó a dudar de la universalidad de este enunciado y se han hecho pruebas donde LOOCV no tiene mayor varianza. Sin embargo sigue siendo computacionalmente más demandante y el beneficio del menor sesgo no era suficiente para darle demasiada importancia.

Cross-Validation en problemas de clasificación.

Los procedimientos vistos hasta ahora son útiles tanto para variables continuas como para problemas de clasificación. Vinimos usando ejemplos donde la medida del error era el MSE (variable dependiente continua) pero podemos aplicar todo de la misma manera utilizando alguna medida de clasificación como la cantidad de observaciones mal clasificadas. Todo el resto se mantiene y es válido, tanto sete de validación, como LOOCV o K-Fold.

Bootstrap

El bootstrap es una herramienta estadística muy extendida que se utiliza para cuantificar la incertidumbre asociada a algún estimador o método de aprendizaje estadístico. Un ejemplo sencilla sería que se puede usar para estimar los errores estándar de los coeficientes de una regresión lineal. Sin embargo lo poderoso de esta herramienta es que es utilizable en muchísimos métodos de aprendizaje, incluso en algunos donde es difícil estimar la varianza o esta no es calculado por los paquetes estadísticos.
Idealmente para estimar la variabilidad de un estimador lo que uno haría es ajustar un modelo n veces y ver cómo varía el estimador a lo largo de esos n modelos utilizando n muestras. Sin embargo no es habitual tener tantos datos ni muestras disponibles. Mismo uno querría utilizar todos los datos en simultáneo posiblemente para reducir el sesgo. Acá es donde bootstrap se luce ya que permite emular el proceso de obtener nuevas muestras de datos a partir de nuestros datos de entrenamiento. En vez de muestrear de manera independiente sobre la población lo que se hace es muestrear n veces con reposición de nuestro set de entrenamiento, generado n muestras a partir de nuestros datos originales. Ya con nuestras nuevas muestras (provenientes todas del dataset original) podemos calcular n modelos y por ende n veces el mismo estimador, pudiendo estimar el desvío estándar de este.
En el fondo lo que se hace es suponer que nuestra muestra es representativa de la población y es nuestra mejor aproximación. Luego obtenemos muestras de estos datos que son nuestra versión reducida de la población. Posiblemente haya algún sesgo pero es una herramienta bastante útil para estimar la variabilidad de nuestros estimadores.

Generamos un ejemplo para ver cómo funciona.

Empezamos generado una población de y que depende x con intercepto 5 y b1 = 5.

library(ggplot2)
set.seed(1)
x <-rnorm(10000, mean = 2, sd = 3)
y <- 4 + 5*x + rnorm(10000,0,4)
df <- cbind.data.frame(y,x)

ggplot(data = df, aes( x =x, y =y )) + 
  geom_point()

Primero vemos el caso ideal que sería poder obtener muchas muestras de la población y ajustar modelos a estas. Luego veremos como varían nuestros coeficientes.

# Muestras de la población
results_pop <- data.frame(b0 = double(), b1 = double())
set.seed(123)
for (i in 1:1000){
  df_train <- df[sample(nrow(df),size = 500,replace = FALSE),]
  ml_train <- lm(formula = y ~ x, data = df_train)
  results_pop[i,1] = ml_train$coefficients[[1]]
  results_pop[i,2] = ml_train$coefficients[[2]]

}

summary(results_pop)

##        b0              b1       
##  Min.   :3.198   Min.   :4.804  
##  1st Qu.:3.834   1st Qu.:4.969  
##  Median :3.974   Median :5.005  
##  Mean   :3.975   Mean   :5.004  
##  3rd Qu.:4.115   3rd Qu.:5.043  
##  Max.   :4.615   Max.   :5.169

print(paste0("El desvío estándar de b0 a partir de 1000 modelos es ", sd(x = results_pop$b0)))

## [1] "El desvío estándar de b0 a partir de 1000 modelos es 0.207882713489026"

ggplot(data = results_pop) + 
  geom_histogram( aes( x = b0), fill = "white", colour = "black")  + 
  geom_vline( aes(xintercept = mean(b0)), colour = "red", size = 1)

Vemos que estimando 1000 modelos a partir de 500 observaciones independientes de la población original obtenemos para b0 estimaciones centradas aproximadamente en el valor real (3.975) pero con un mínimo encontrado en 3.198 y un máximo en 4.615. El desvío estándar de la estimación es de 0.2078. A su vez mostramos un histograma de cómo se distribuye la estimación de b0.

Ahora simulemos un caso real donde solo tenemos una muestra de 500 observaciones y es todo con lo que podemos trabajar. Como primera medida estimamos una regresión lineal y vemos qué parámetros ajustan mejor nuestros datos.

# Muestras de la población
results_sample <- data.frame(b0 = double(), b1 = double())
set.seed(456)
df_train_sample <- df[sample(nrow(df),size = 500,replace = FALSE),]
ml_train_sample <- lm(formula = y ~ x, data = df_train_sample)

results_sample[1,1] <- ml_train_sample$coefficients[[1]]
results_sample[1,2] <- ml_train_sample$coefficients[[2]]

knitr::kable(results_sample, caption = "Coefficients")

Table: Coefficients

b0	b1
3.89621	5.025914

Vemos que a partir de entrenar el modelo con las 500 observaciones obtenemos un intercepto de 3.896 y un b1 estimado de 5.026. Nosotros, como conocemos la población, sabemos que el intercepto no es del todo preciso ya que el real es 4 sin embargo en un caso real eso no lo sabríamos. Nos interesaría saber qué variabilidad tiene ese coeficiente para tener una medida de qué tan variable es nuestro resultado.
Para una regresión lineal eso se puede saber ya que no es difícil calcular la varianza de los estimadores, pero con modelos más complicados no siempre se puede y ahí es donde bootstrap ayuda realmente. Acá lo hacemos con la regresión lineal porque es lo más sencillo de mostrar.
Suponiendo que queremos obtener una estimaación de la variabilidad del coeficiente estimado b0 procedemos con bootstrap.

Fijense que lo que hacemos es distinto al primer caso. Acá tomamos 10000 muestras no de la población sino de nuestro set de 500 observaciones. Estas muestras son también de 500 observaciones, la diferencia es que es con reposición por lo tanto una misma observación puede figurar más de una vez.

# Muestras de la población
results_bootstrap <- data.frame(b0 = double(), b1 = double())
set.seed(789)
for (i in 1:10000){
  df_train_bs <- df[sample(nrow(df_train_sample),size = 500,replace = TRUE),]
  ml_train_bs <- lm(formula = y ~ x, data = df_train_bs)
  results_bootstrap[i,1] = ml_train_bs$coefficients[[1]]
  results_bootstrap[i,2] = ml_train_bs$coefficients[[2]]

}

summary(results_bootstrap)

##        b0              b1       
##  Min.   :3.245   Min.   :4.791  
##  1st Qu.:3.822   1st Qu.:4.982  
##  Median :3.962   Median :5.020  
##  Mean   :3.962   Mean   :5.019  
##  3rd Qu.:4.100   3rd Qu.:5.057  
##  Max.   :4.709   Max.   :5.205

print(paste0("El desvío estándar de b0 a partir de 10000 modelos es ", sd(x = results_bootstrap$b0)))

## [1] "El desvío estándar de b0 a partir de 10000 modelos es 0.207126704818891"

ggplot(data = results_bootstrap) + 
  geom_histogram( aes( x = b0),fill = "white", colour = "black") + 
  geom_vline( aes(xintercept = mean(b0)), colour = "blue", size = 1)

Voilà. Corrimos 10000 iteraciones de nuestro modelo a partir de 10000 muestras de nuestra data original. El desvío estándar de b0 para bootstrap quedó de 0.2071. Que si comparamos con el de 1000 muestras independientes que era 0.2078 es prácticamente igual. A su vez, podemos calcular el desvío teórico de b0 a partir del modelo (la solución fácil).

summary(ml_train_sample)

## 
## Call:
## lm(formula = y ~ x, data = df_train_sample)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.2193  -2.8265   0.0281   2.7421  10.7577 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.89621    0.21607   18.03   <2e-16 ***
## x            5.02591    0.05902   85.16   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.98 on 498 degrees of freedom
## Multiple R-squared:  0.9357, Adjusted R-squared:  0.9356 
## F-statistic:  7252 on 1 and 498 DF,  p-value: < 2.2e-16

Vemos que desde R el modelo nos devuelve que b0 tiene un desvío de 0.21607. Prácticamente igual al desvío de las muestras independientes como al de bootstrap. Por otra parte vemos que el promedio de b0 estimado en bootstrap es mucho más cercano a 4 que el estimado con una sola iteración y quedó mucho más cerca que el promedio de los estimados mediante muestras independientes. Nada mal no?

2019-05-25
in estadistica, R
13 min read

Regresion Lineal - ISLR Capítulo 3

La regresión lineal simple es un método muy directo para estimar una variable cuantitativa Y en base a un solo predictor X. Asume que hay una relación lineal entre X e Y. $$ Y \approx \beta_0 + \beta_1X$$ $\beta_0$ y $\beta_1$ son dos constantes desconocidas que representan al intercepto y a la pendiente del modelo lineal. Son los coeficientes o parámetros. Con nuestros datos podemos estimar coeficientes para predecir futuros valores de Y basados en X y nuestro modelo.

Estimación de Coeficientes

Los coeficientes que buscamos son $\hat \beta_0$ y $\hat \beta_1$ (estimados, por eso el sombrero) son aquellos que generen una recta que pase lo más cerca posible de todos nuestros datos de entrenamiento. Hay varias manera de definir "cerca" pero la más usada es el enfoque de mínimos cuadrados.

Supongamos un ejemplo donde tenemos datos de horas trabajadas por ciertos individuos y la paga que reciben. Supongamos a fines del ejemplo que la relación entre salario y horas es lineal (sabemos que no es real...)

Cuando estimemos $\hat \beta_0$ y $\hat \beta_1$ obtendremos después un valor $\hat y_i$ para cada valor de $x_i$ (cada observación), que será el resultado de la predicción de nuestro modelo para ese valor de horas trabajadas. $$ \hat y_i = \hat \beta_0 + \hat \beta_1 x_i$$ Luego $e_i = y_i - \hat y_i$ representa el residuo, que es la diferencia entre el valor real del salario para esa observación y el valor que predice nuestro modelo.
Una métrica importante a saber es la suma de resiudos al cuadrado (RSS por siglas en inglés) que es: $$ RSS = e_1^2 + e_2^2 + ... + e_n^2$$ o de manera equivalente: $$ RSS = (y_1 - \hat \beta_0 - \hat \beta_1x_1) ^2 + (y_2 - \hat \beta_0 - \hat \beta_1x_2) ^2 + ... + (y_n - \hat \beta_0 - \hat \beta_1x_n) ^2 $$

Que basicamente es la suma de todas las diferencias entre lo predicho por nuestro modelo y el dato real de nuestro set, elevadas al cuadrado. Esto último es principalmente para evitar que se compensen los errores. Sobreestimar por 10 y luego subestimar por 10 tiene como suma de errores 0. Si elevamos esas diferencias al cuadrado, todos los errores serán positivos y se acumularán. En este caso seria $10^2$ + $(-10)^2$, que es 200.
El enfoque de mínimos cuadrados estima $\hat \beta_0$ y $\hat \beta_1$ de tal manera que el RSS sea el mínimo posible dados los datos.

Usando un poco de cálculo se puede demostrar que los parámetros que minimizan RSS son: $$ \hat \beta_1 = \frac{\sum_{i = 1}^n (x_i - \bar{x}) (y_i - \bar{y})}{\sum_{i=1}^n (x_i - \bar{x})^2}$$ $$ \hat \beta_0 = \bar{y} - \hat \beta_1\bar{x}$$ donde $\bar{y}$ y $\bar{x}$ son las respectivas medias muestrales.

En nuesto caso usando este set de datos generado ficticiamente obtenemos $\hat \beta_0$ = 2.000713 × 10⁴
y $\hat \beta_1$ = 299.82

Recordemos que esta es una estimación en base a los datos y no sabemos los verdaderos parámetros de la DGP real(proceso generador de datos). En este caso yo si lo sé porque generé los datos pero en la vida real es inaccesible.
Lo que hicimos fue estimar, a partir de un set de datos, ciertos coeficientes o característica de una población mucho más amplia. (Todos los trabajadores del país..)

Precisión de nuestros estimadores

En este caso teníamos una sola muestra pero podríamos haber tenido muchas muestras (K) de la misma población (muchos sets de datos con horas trabajadas y salarios).
Si estimáramos los coeficientes para cada uno de esos sets obtendríamos K pares de coeficientes, cada uno calculado con las particularidades de esos sets.
Se puede demostrar que el promedio de una cantidad grande de estimadores provenientes de muchas muestras se centra en el verdadero valor poblacional (si el modelo es correcto). Es decir que el promedio de los K $\hat \beta_1$ va a centrarse en el verdadero valor poblacional de $\beta_1$ ( y lo mismo para $\beta_0$).
Pero estos K parámetros centrados en el verdadero valor van a tener cierta dispoersión, es decir, pueden estar todos muy cerca del verdadero o estar muy dispersos pero que en promedio si quede centrado. Esto determina que tan preciso es el coeficiente que estimemos de una muestra. Este desvío estándar de los parámetros (SE) puede estimarse y depende de la varianza del error del modelo.

Puede ser útil para calcular los intervalos de confianza de los parámetros. Estos son intervalos que con X% de probabilidad contienen al verdadero valor del parámetro poblacional. Lo más habitual es calcular el intervalo de confianza al 95%. Para $\hat \beta_1$ esto es aproximadamente: $$ \hat \beta_1 \pm 2 \cdot SE(\hat \beta_1)$$ La interpretación sería que de 100 intervalos que construya de esta manera (de 100 muestras distintas), 95 van a tener al verdadero valor de $\beta_1$.

Por otra parte podemos realizar un test de hipotésis de los coeficientes. El más común es testear la siguiente hipótesis nula:
H_0 : No hay relación entre X e Y
contra la hipótesis alternativa
H_1 : Hay alguna relación entre X e Y
Lo cual se traduce en: $$ H_0 : \beta_1 = 0 $$ $$ H_1 : \beta_1 \neq 0$$

Lo que se hace es determinar si $\hat \beta_1$ está lo suficientemente lejos de 0 como para rechazar la hipótesis nula. Qué tan lejos es suficiente depende en gran parte del desvío estándar (SE) del coeficiente. Si el SE es grande , necesitaremos valores elevado de $\hat \beta_1$ para estar tranquilos con que el valor real no puede ser 0.
Para esto lo que se hace es calcular el estadístico t: $$ t = \frac{\hat \beta_1 - 0}{SE(\hat \beta_1)}$$ Que mide cuantos desvíos estándar $\hat \beta_1$ está alejado de 0. Si no hay relación entre X e Y se espera que el estadístico tenga una distribución t con n - 2 grados de libertad. Dado ese supuesto, lo que se hace es calcular la probabilidad de obtener un valor de t como el de nuestro estadístico, si este proviene de una distribución t con n-2 grados de libertad. Esta probabilidad se la conoce como p valor. Sería qué tan probable es encontrar un valor al menos tan grande como el de t si este proviniera de la distribución t con n-2 grados de libertad. Si esta probabilidad es muy chica (el umbral habitual es 0.05 pero depende del trabajo) uno rechaza la hipótesis nula en favor de la alternativa, suponiendo que sí hay una relación entre X e Y.

Precisión del modelo.

Naturalmente a uno le interesa saber qué tan bien ajusta nuestro modelo a los datos.

El método más habitual para regresión lineal es el R^2. Toma valores entre 0 y 1 porque es la proporción de la varianza de Y explicada por nuestro modelo.

\[ R^2 = \frac{TSS - RSS}{TSS} = 1 - \frac{RSS}{TSS}\]

donde $TSS = \sum (y_i -\bar{y})^2 $ es la suma de cuadrados totales y RSS es la suma de errores cuadrados que ya definimos antes. TSS mide la varianza total de Y y representa la variabilidad total inherent de la variable dependiente antes de correr la regresión. Por el contrario RSS mide la variabilidad que queda sin explicar por nuestro modelo (recuerden que proviene de los residuos). Por lo tanto el numerador TSS - RSS mide la parte de la variabilidad de Y que sí pudo ser explicada por el modelo, y lo divide por la variabilidad total. $R^2$ mide entonces la proporción de la variablidad que pudo ser explicada usando X. Cuanto más cerca de 1, mejor.

Regresión con Múltiples Predictores

Suena mucho más lógico tratar de explicar una variable dependiente no solo por una independiente si no por varias. La regresión lineal simple puede ampliarse a regresión lineal múltiple donde nuestro modelo pasa a ser: $$ \hat y_i = \hat \beta_0 + \hat \beta_1 x_i + ... + \hat \beta_p x_p$$ Y mantenemos un término de error con distribución normal y media 0.

En esencia la idea es la misma, explicar la variabilidad de Y basado en la variabilidad de nuestros predictores. La metodología para estimar los coeficientes suele ser Mínimos Cuadrados como vimos antes, sin embargo la solución no suele ser tan fácil de expresar y es más sencillo verlo en términos matriciales o simplemente ver los resultados desde el programa estadístico que estemos usando. No olvidar que varias regresiones simples no pueden sumarse en una resgresión múltiple, es decir, los coeficientes de las regresiones simples no tienen por que ser los mismos ni por qué mantener el signo cuando se juntan todas las variables en un solo modelo. Esto sucede porque la regresión múltiple estima coeficientes controlando por todas las otras variables, es decir, quitando el efecto de las otras. Por eso es que por separado quizás dos variables son significativas pero en una regresión múltiple solo una de ellas lo es. En general esto viene dado porque están correlacionadas y se comportan de manera similar. Al final del día no es fácilmente distinguible cuál es realmente la que lidera el efecto.

Hay relación entre la dependiente y los predictores? (test F)

En regresión simple vimos el test de hipótesis para ver si el coeficiente de $\hat \beta_1$ era significativamente distinto de 0. En regresión múltiple lo que debemos hacer es chequear si todos nuestros coeficientes son distintos de 0. (Y no uno por uno)

Lo cual se traduce en: $$ H_0 : \beta_1 = \beta_2 = ... = \beta_p = 0 $$ $$ H_1 : \text{Al menos algún } \beta_j \neq 0$$ El test de hipótesis se hace calculando el estadístico F. $$ F = \frac{(TSS - RSS) / p}{RSS / (n - p -1)}$$

Si los supuestos del modelo lineal se cumplen puede probarse que $E[RSS/(n - p-1)] = \sigma^2$ y si $H_0$ es verdadera $E[(TSS-RSS)/p] = \sigma^2$. Por lo tanto si no hay relación el estadístico F será cercano a 1 y si en realidad la hipótesis alternativa es verdadera el numerador será mayor que $\sigma^2$ y por lo tanto F será mayor que 1. Dependiendo de n, p y del nivel de significatividad que busquemos F deberá superar un umbral distinto para poder rechazar la hipótesis nula.
Es inevitable mirar los p-valores individuales sin embargo debemos tener cuidado particularmente cuando tenemos muchas variables. Por definición, algunos coeficientes saldrán significativos por azar aunque no tengan relación con la variable dependiente. En el caso típico de significativdad del 95%, esto sucede en promedio el 5% de las veces. Con muchas variables nuestras posibilidades de encontrarnos con al menos algún falso significativo aumentan notoriamente por lo que hay que mirar con cuidado. Por su parte, el estadístico F corrige en su cálculo por la cantidad de coeficientes y por lo tanto no se ve afectado por este problema.

Selección de variables

Cuando tenemos un set de datos grande es habitual tener que seleccionar cuáles son las variables importantes para el modelo. Más allá del conocimiento del problema (fundamental), idealmente lo mejor es probar una gran cantidad de modelos y con alguna métrica de comparación seleccionar los mejores. El problema es que la cantidad de modelos posible crece exponencialmente con la cantidad de variables y esto no es posible.
En el libro los autores mencionan como alternativas Forward Selection, Backward Selection y Selección mixta. Básicamente son enfoques que prueban una muestra de todos los modelos posibles según la significatividad de las variables. Son métodos iterativos. Habiendo avanzado la disciplina, llegado el caso buscaría otros métodos vigentes para atacar este problema.

Ajsute del modelo

Para verificar el ajuste del modelo se sigue usando el $R^2$ como métrica principal. En este caso es equivalente a la correlación al cuadrado de Y e $\hat Y$. Un punto a tener en cuenta es que el R^2 nunca puede disminuir al agregar variables ya que el peor escenario posible es que la nueva variable tenga coeficiente de 0 y el ajuste quede igual que antes. Lo que se hace para controlar por esto y poder comparar modelos es ajustar el $R^2$ por la cantidad de variables utilizadas o usar el RSE. De cualquier manera lo importante es recordar que el R2 sigue siendo útil en la regresión lineal múltiple.

Otras consideraciones

Las variables independientes admiten variables categóricas! (Binarias o multiclase). Ej: Educación máxima alcanzada. Hay que mirar con atención la interpretación. Alteran el intercepto según la categoría de la observación y puede alterar pendientes si se las incluye en interacción con alguna variable continua.
El modelo que venimos viendo el aditivo y lineal, pero podemos remover esos supuestos. Por ejemplo podemos agregar interacción entre variable y por lo tanto relajar la aditividad. Esto significa que las variables se modelan multiplicadas entre sí por ejemplo.

\[ Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_1X_2 + \epsilon \]

Podemos aproximar relaciones no lineales extendiendo el modelo a regresión polinómica. Suponiendo que los datos provienen de un modelo polinómico, podemos ver en el siguiente gráfico cómo cambia al agregar el término de no linealidad. La línea naranja es la regresión lineal simple y la celeste (que ajusta casi perfecto) es la regresión polinómica que respeta el proceso generador de los datos (la ecuación que se ve en el gráfico). Vemos que la variable Y depende de X linealmente pero también de X al cuadrado, lo que le da la curvatura.

Potenciales problemas

No linealidad de los datos

Si la relación entre nuestras variables independientes y la dependiente no es lineal nuestro modelo va a tener sesgo alto. Para el modelo simple es fácil de ver al graficar X vs Y en un gráfico de puntos pero con muchas variables eso ya no es tan sencillo. Un buen enfoque es realizar una regresión lineal y graficar los residuos contra los valores predichos.

Es un caso un poco extremo pero supongamos que la relación es no lineal, polinómica de orden 2 (como el ejemplo de arriba).
Si nosotros corremos una regresión lineal, nuestros residuos van a seguir un patrón muy obvio.

Claramente en ese gráfico el residuo no está centrado en 0... En cambio, si nosotros corremos una regresión lineal para un modelo realmente lineal, o para este caso, corremos un modelo no lineal, deberíamos ver una nube de puntos dispersa para los resiudos, centrada en 0 y con algún desvío estándar. Idealmente veríamos algo asi.

Correlación de los términos de error

Uno de los supuestos de la regresión lineal es que los errores no están correlacionados, es decir que el error $\epsilon_i$ de una observación no nos aporta información acerca del error $\epsilon_j$ de otra observación. Son independientes.
Si esto no se cumple lo que sucede es que el SE de los coeficientes estimados es menor al real y puede llevarnos a confiar más en nuestro modelo de lo que deberíamos.
Los errores correlacionados suelen suceder más frecuentemente en series de tiempo pero también pueden darse en estudios experimentales mal diseñados.

Varianza de los términos de error no constante. (Heterocedasticidad)

Otro supuesto de la regresión lineal es que la varianza de los errores es constante $Var(\epsilon_i) = \sigma^2$. Esto no siempre es el caso. En este ejemplo vemos como los resiudos siguen centrados en 0 pero con una dispersión mucho mayor a medida que avanzamos en el eje X. Entre las soluciones para este problema se encuentra transformar la variable dependiente - $ln(y)$ por ejemplo, o utilizar Mínimos cuadrados ponderados, que pondera por la inversa de la varianza. El libro no se explaya mucho más al respecto en este capítulo.

Outliers

Los outliers son observaciones cuya variable dependiente tienen valores que se alejan mucho del patrón regular de los datos, por ejemplo debido a un error de medición o problema al registrar la información. Los outliers pueden tener diversas consecuencias en los modelos lineales. Puede afectar la estimación de los parámetros, puede afectar el ajuste del modelo (caída del $R^2$) o puede por ejemplo aumentar los intervalos de confianza ya que el outlier afecta el RSE que es común a todos los intervalos. Todo esto puede ser generado por una sola observación. Generalmente si no se distinguen en el análisis exploratorio pueden saltar a la vista analizando los resiudos del modelo (o los residuos estandarizados).

Puntos con alto "leverage"

Si alguno tiene una traducción satisfactoria bienvenido sea.. Los puntos con alto leverage son aquellos cuyo valor de la variable independiente se aleja del rango estándar. Las observaciones con esta característica tienden a afectar en buena medida a la curva ajustada y por ende a los parámetros de nuestro modelo. Nuestra estimación por mínimos cuadrados puede verse muy influenciada por estos puntos e invalidar el ajuste por eso es muy necesario identificar estas observaciones.
En regresión simple es sencillo de ver porque resaltan si graficamos una nube de puntos pero en regresión múltiple es más difícil de ver ya que debemos encontrar anomalías en el conjunto de todas las variables. Es decir que una observación puede estar en el rango individual de cada variable pero si miramos a nivel conjunto, esa combinacion dentro de los rangos individuales es súper anómala. Con más de dos variables independientes se dificulta identificar visualmente. Para ayudar en estos casos se puede calcular el estádistico de leverage en algún programa estadístico.

Colinealidad

Este problema refiere a la alta correlación entre variables independientes del modelo, es decir que tienden a aumentar o decrecer de manera conjunta. Esto genera que sea difícil (o imposible en el extremo) diferenciar el impacto de cada una de ellas en la variable dependiente.
En una regresión lineal esto se traduce en aumento de la varianza de los estimadores y por ende incertidumbre sobre los parámetros estimados. A modo intuitivo, con variables con alta correlación puede haber una gran cantidad de combinaciones de coeficientes para estas variables que resulten en un mismo ajuste ($R^2$) y por ende mínimos cuadrados es indistinto frente a ellos. Cambiando alguna observación puede que el modelo pase de una combinación a otra muy disinta en ese arco de posibilidades. Otra consecuencia es que el aumento de la varianza de los coeficientes reduce el estadístico t que miramos para la significatividad y puede que lleve a no rechazar una hipótesis nula que debía ser rechazada. La potencia del test de hipótesis se ve disminuida por la colinealidad.
No solo la correlación sirve para detectar colinealidad. Puede existir multicolinealidad en donde varias variables son colineales aún sin tener alta correlación de a pares. Posiblemente se deba a combinación lineal generada por algunas de las variables. Para estos casos lo que se puede mirar es el VIF ( Variance Inflation Factor) en inglés. Este estadístico se calcula para cara variable y compara la varianza del estimador al tener la variable en el modelo versus ajustando un modelo solo con esa variable. Cuanto mayor es el VIF, mayores problemas de colinealidad resalta. Se puede calcular con la siguiente formula donde $R^2_{X_j|X_{-j}}$ es el $R^2$ de la regresión de $X_j$ contra todas las otras variables independientes del modelo. $$ VIF(\hat \beta_j) = \frac{1}{1 - R^2_{X_j|X_{-j}}}$$ La solución a este problema suele ser descartar alguna de las variables o agruparlas de alguna manera para quedarnos con una única variable que represente a ambas.

2019-05-06
in estadistica
8 min read

Aprendizaje Estadístico - ISLR Capitulo 2

Suponemos que las variables que encontramos en un set de datos son generadas a través de un proceso generador de datos (DGP por sus siglas en inglés) cuya expresión es: $$ Y = f(X) + \epsilon $$

Donde Y es la variable, en este caso la dependiente o la que queremos explicar. f(X) es una función respecto a otra/s variable/s X (independientes) y $\epsilon$ es el error irreducible, es decir un valor aleatorio con media 0 pero que no depende de otras variables, es al azar. Puede referir a errores de medición, cambios inmesurables en las situaciones del experimento o simplemente azar en la generación real de los datos. Cabe destacar que f(X) es desconocida para nosotros y justamente lo que queremos explorar con el análisis estadístico. Puede tenerse suposiciones o conocimiento de la forma funcional (lineal, no lineal, etc) pero en principio no tenemos mayores certezas y esperamos aprender a partir de la muestra que analizamos.

Por qué estimar f(X)?

Los dos principales motivos para interesarse en f(X) son predicción (de Y) e inferencia de los parámetros de f(X).

Predicción

Queremos predecir valores de Y para nuevos datos X. Como $\epsilon$ en promedio es 0 podemos aproximar Y de la forma: $$ \hat Y = \hat f(X)$$ La precisión de $\hat Y$ va a depender del error reducible y del error irreducible. El primero depende de qué tan bien nos aproximemos a la verdadera f(X) y puede ser potencialmente reducido si utilizamos las técnicas más adecuadas para el caso. El segundo error es justamente irreducible y es porque nuestra aproximación no puede tener en cuenta a $\epsilon$. El término aleatorio introducido por esa variable no lo podemos estimar para cada observación y por lo tanto debemos convivir con ese margen de error.

Suponiendo que tenemos una estimación $\hat f$ y un set de datos X puede probarse que: $$ E(Y - \hat Y)^2 = E[f(X) + \epsilon - \hat f(X)]^2$$ $$ = \underbrace{[f(X) - \hat f(X)]^2}\text{Reducible} + \underbrace{Var(\epsilon)}\text{Irreducible}$$

Donde $E(Y-\hat Y)^2$ es el promedio o valor esperado de la diferencia al cuadrado del valor real de Y y de la predicción correspondiente. $Var(\epsilon)$ es la varianza del término de error $\epsilon$.

Inferencia

Este enfoque se basa en entender la relación entre las variables de X y la dependiente. Es necesario entender bien la f(X) elegida para poder interpretar sus coeficientes y poder ver qué variables están asociadas con Y, cómo es esa relación, cuál es la forma de la función f(X), etc para poder actuar sobre las variables X o comprender su efecto aunque no siendo tan exigentes con el poder de predicción de nuestro modelo.

Como Estimar f(X)?

Métodos Paramétricos

Los métodos paramétricos se conforman por dos etapas.
La primera es asumir o suponer la forma funcional de f(X). Podes definir por ejemplo que f(X) es una función lineal de la forma $$ f(X) = \beta_0 + \beta_1X_1 + \beta_2X_2 + ... + \beta_pX_p$$

Una vez definida la forma del modelo la segunda etapa consiste en estimar los parámetros con algún método, a partir de los datos de entrenamiento. En este caso sería estimar todos los $\beta$. Por ejemplo para las funciones lineales se suele utilizar el método de mínimos cuadrados ordinario.
La ventaja de definir una forma funcional es que luego es más sencillo estimar sus parámetro y el problema se reduce a eso finalmente. Por el otro lado, posiblemente la forma que elijamos no sea exactamente igual a la real (DGP) y tengamos que aceptar que va a haber errores debido a eso. Si estamos muy lejos de la forma real esos errores serán groseros. Existen modelos flexibles que permiten ajustar modelos con diferentes formas de f(X) pero en general requieren estimar más parámetros y son más propensos a sufrir sobreajuste/overfitting que básicamente es ajustarse mucho al ruido o error ($\epsilon$) en los datos de entrenamiento y luego ajustar mal en testeo.

Métodos No Paramétricos

Los métodos no paramétricos no requieren definir explícitamente una forma funcional de f. Buscan un f que sea lo más cercano posible a los datos sin ser demasiado estricto o flexible. Al no asumir una forma puede cubrir potencialmente un rango mucho mayor. El problema es que al no reducir el problema a estimar parámetros necesitan muchas más observaciones para estimar f de forma medianamente precisa. En general uno tiene que decidir el nivel de "suavidad" del modelo, lo cual afecta que tan variable termina siendo la estimación. Sirve para encontrar el punto de fleixibilidad/rigidez del modelo que queremos para que no sobreajuste (ni falle demasiado).

Predicción vs Interpretabilidad

Uno puede elegir entre modelos flexibles o más rígidos. Con pocas observaciones a veces uno no puede alejarse mucho de los rígidos pero suponiendo que uno tiene muchos datos, a veces puede igualmente elegir rígidez frente a modelos flexibles que permiten ajustar varias formas de f. El motivo es que generalmente los modelos restrictivos son más fáciles de interpretar y se le puede dar un significado claro a sus coeficientes mientras que con formas muy flexibles no es sencillo entender el impacto de las variables de manera individual. La elección va a depender del objetivo del análisis y de qué tan bien o mal nuestros modelos ajustan a los datos.

Modelos Supervisados vs No Supervisados

Nuestros datos pueden tener una variable dependiente que queremos explicar o predecir en base a un set de variables independientes, con algún modelo a definir. Los casos de este estilo son llamados supervisados porque sabemos la "respuesta" (nuestra variable Y) y podemos validar nuestros modelos contra la realidad.
Si los datos no tienen una variable dependiente lo que se puede hacer es un análisis no supervisado donde por ejemplo lo que se puede hacer es agrupar las observaciones en clusters o grupos. Es decir segmentar en distintas clasificaciones y descubrir patrones. El desafío es que no hay en los datos nada contra qué validarlo, aunque sí contra el conocimiento del dominio o de la temática.

Regresión vs Clasificación

En los modelos supervisados, nuestra variable dependiente puede ser cuantitativa o cualitativa.
En el primer caso la variable toma valores númericos, como por ejemplo la altura de una persona, el precio de una propiedad, etc. Son problemas de regresión.
En el segundo caso la variable dependiente puede tomar el valor de una clase o categoría. Por ejemplo, género de una persona, si paga o no paga su deuda, etc. Son problemas de clasificación.

Midiendo el ajuste del modelo

Regresión

Para problemas de regresión una de las medidas más utilizadas es el Error Cuadrático Medio (MSE por sus siglas en inglés). $$ MSE = \frac{1}{n} \sum_{i=1}^n(y_i - \hat f(x_i))^2 $$

Es básicamente la diferencia promedio entre la realidad y lo que predice nuestro modelo elevado al cuadrado. Esto último es para que los errores sean siempre positivos aunque subestimemos o sobreestimemos (y por su comodidad para cálculos matemáticos).

En primero lugar se calcula este valor con los datos de entrenamiento sin embargo lo que realmente importa es como performa el modelo en datos de testeo, es decir en datos que no fueron utilizados para estimar f(X). Podemos decir que cada modelo debería tener un MSE de entrenamiento y un MSE de testeo. Debido a la posibilidad de sobreajuste y a las diferencias en muestras nada garantiza que el modelo que estimemos con menor MSE en entrenamiento también sea el de menor MSE en testeo.

A medida que aumentamos la flexibilidad de un modelo (sus grados de libertad) el MSE en entrenamiento va a disminuir, ya que tiene más herramientas para ajustarse a los datos pero puede que sobreajuste y por lo tanto no se traduzca en un menor MSE en testeo.

El tradeoff entre Sesgo y Varianza

No está demostrado en el libro pero es posible descomponer el MSE esperado de una observación de testeo en sesgo de $\hat f(x_0)$, varianza de $\hat f(x_0)$ y varianza del error irreducible $\epsilon$.

\[ E(y_0 - \hat f(x_0))^2 = Var(\hat f(x_0)) + [Sesgo(\hat f(x_0))]^2 + Var(\epsilon)\]

El lado izquierdo de la ecuación es el MSE esperado y coresponde MSE de testeo promedio que obtendríamos si estimaramos f utilizando una gran cantidad de sets de entrenamiento y testearamos cada uno en $x_0$.

Algunas observaciones:

El MSE nunca puede ser menor que la varianza de $\epsilon$. Es un término fijo y por eso se lo llama error irreducible.
La varianza es cuanto cambiar $\hat f$ si utilizamos otro set de entrenamiento. Siempre va a cambiar con otro set pero idealmente ese cambio no debería ser grande. Modelos muy flexibles tienden a cambiar más frente a distintos sets y son más inestables.
Sesgo es el error provocado por la diferencia entre el modelo elegido y el verdadero proceso generador de los datos. En general modelos más flexibles tienen menor sesgo ya que pueden ajustar mayor variedad de formas funcionales.
Al aumentar la flexibilidad de un modelo en general reducimos el sesgo pero aumentamos la varianza. En general en un primer momento el sesgo suele disminuir a mayor velocidad de lo que aumenta la varianza y por lo tanto el MSE esperado se reduce. Sin embargo llega un punto donde mayor flexibilidad reduce menos el sesgo que lo que aumenta la varianza y el MSE empieza a aumentar. Es el primer indicio de sobreajuste. Por eso se habla de tradeoff o "balance".
En la realidad donde la verdadera f del DGP es inobservable no suele ser posible calcular explícitamente el MSE de testeo, el sesgo o la varianza de un método estadístico pero el proceso de fondo aplica y siempre debemos tener en mente el tradeoff.

Clasificación

Para problemas de clasificación uno de los enfoques más frecuentes para cuantificar la precisión de una función estimada $\hat f$ se suele usar el porcentaje de error en los datos de entrenamiento.

$$ \frac{1}{n} \sum_{i=1}^nI(y_i \neq \hat y_i) $$ Básicamente es el porcentaje de observaciones clasificadas erroneamente. Al igual que con MSE es de gran importance el porcentaje de error en los datos de testeo.

Clasificador de Bayes

No lo demuestra en el libro pero la mejor manera de reducir el porcentaje de error en test es asignar a cada observación la clase con mayor probabilidad (según el DGP ). Es un concepto muy sencillo, dado X, asignar la clase cuya chance de acierto sea mayor.

$$ Pr(Y = j | X = x_0) $$ El porcentaje de error de Bayes (es decir el error luego de clasificar siguiendo esa regla) es análogo al Error Irreducible de regresión. Hay que tener en cuenta que la distribución condicional de Y dado X no lo sabemos en los casos aplicados en la vida real, sería como saber la función f(X) o el DGP y por lo tanto no lo podemos calcular.

K vecinos más cercanos (KNN en inglés)

Idealmente uno querría aplicar el clasificador de Bayes pero es imposible ya que no sabemos la distribución real de los datos (es justamente lo que queremos estimar). KNN intenta aproximarse a la distribución condicional para clasificar las observaciones. Lo que hace este método es, dado un valor de K que elegimos nosotros, clasificar cada nueva observación según la clase mayoritaria entre las K observaciones más cercanas a esta.

$$ Pr(Y = j | X = x_0) = \frac{1}{K} \sum_{i \in N_0} I(y_i = j) $$ El valor que seleccionemos de K afecta en gran medida las predicciones del modelo. Un K menor hace más variable el modelo ya que selecciona menos observaciones y por lo tanto pocos cambios en el set de entrenamiento cambian la clasificación. Suele reducir el sesgo pero ser mas variable. Es análogo a hacer más flexible un modelo en regresión. Valores de K más grandes seleccionan puntos en un entorno más abarcativo y por lo tanto suele ser más constante pero con sesgo superior.
Al igual que en regresión hay que tener cuidado con el sobreajuste. Reducir K garantiza menos errores en los datos de entrenamiento pero pasado un umbral la varianza aumenta en mayor medida y el porcentaje de error en test se incrementa.

Conclusión: Tanto en Clasificación como en Regresión la elección del nivel de flexibilidad es central en el éxito de método de aprendizaje estadístico.

2019-04-03
in estadistica, R
4 min read

Funciones de Probabilidad y Distribucion

Variables Aleatorias

Consideremos un experimento cuyo espacio muestral denominaremos S.
Una funcion valuada en el dominio de los reales definida en S es una variable aleatoria.

En otras palabras es una función que asigna a cada resultado posible de un experimento un valor real.

Por ejemplo:

Si el experimento es lanzar una moneda 10 veces hay 2¹⁰ combinaciones posibles de caras (o) y cruz (x).
Si definimos la variable aleatoria X como cantidad de caras entonces X(s) será la cantidad de caras del experimento.
Si s resulta ser la secuencia ooxxxoxxxo entonces X(s) = 4.

Distribucion de una variable aleatoria

Si tenemos la distribución de probabilidad del espacio muestral del experimento podemos determinar la distribución de probabilidad de cualquier variable aleatoria válida.

Volviendo al ejemplo de la moneda. Dijimos que hay 2¹⁰ combinaciones de cara o cruz. La cantidad de combinaciones de X caras en 10 lanzamientos es $P(X = x) = \binom{n}{x} \frac{1}{2^{10}}$ para $x = 0,1,2,..,10$

Distribuciones Discretas

Una variable aleatoria tiene una distribución discreta si solo puede tomar valores de una secuencia (generalmente finita pero puede no serlo).

La función de probabilidad le otorga una probabilidad puntual a cada valor de esa secuencia.
Los valores por fuera de la secuencia tienen probabilidad = 0
La suma de todas las probabilidades tiene que ser 1

Distribución Uniforme

En el caso de la dsitribución uniforme, supongamos que la variable puede tomar valores de 1 a k. La función de probabilidad será $f(x) = \frac{1}{k}$ para x = 1,2,...,k. Y 0 para todos los otros valores.

si k = 10
Los valores de la variable serán cualquier entero entre 1 y 10
Cada valor tendrá probabilidad $\frac{1}{10}$

Distribución Binomial

En el caso de la dsitribución binomial se asumen dos posibles resultados, uno con probabilidad p y su contraparte con probabilidad 1-p.
Por ejemplo la probabilidad p de que una máquina genere un producto defectuoso y 1-p de que sea no defectuoso.
Si una máquina produce n productos va a generar X productos defectuosos. La variable aleatoria X tendrá una distribución discreta y sus posibles valores irán de 0 a n.
Para cualquier valor de x (entre 0 y n), la probabilidad de que la máquina genere x productos defectuosos entre los n producidos (de una secuencia particular) es $p^{x}q^{(n-x)}$
Como existen $\binom{n}{x}$ distintas secuencias posibles con x defectuosos entre los n productos tenemos que:
$Pr(X = x) = \binom{n}{x}p^{x}q^{(n-x)}$
La función de probabilidad será $f(x) = \binom{n}{x}p^{x}q^{(n-x)}$ para x = 0,1,2,...,n. Y 0 para todos los otros valores.

Para usar esta distribución en R tenemos los siguientes comandos:

Para generar n escenarios al azar donde se producen size productos con probabilidad p de ser defectuosos. El resultado es la variable x por escenario. Es decir la cantidad de defectuosos.
En el primer escenario x = 0, en el segundo x = 1 y así.

set.seed(1)
rbinom(n = 10, size = 5, p = 0.2 )

##  [1] 0 1 1 2 0 2 3 1 1 0

# random binomial

Para saber la probabilidad de obtener x productos defectuosos si una máquina produce size productos y la probabilidad de que produzca un defectuoso es prob.
Hay probabilidad de 0.0264 de obtener 5 defectuosos si producimos 10 con probabilidad 0.2.

dbinom(x = 5, size = 10, prob = 0.2)

## [1] 0.02642412

Para saber la probabilidad acumulada de obtener q o menos productos defectuosos si la máquina fabrica size objetos, con probabilidad de defecto prob. Hay probabiliad de 0.879 de obtener 3 o menos defectuosos si la máquina produce 10 objetos con probabilidad 0.2 de defecto. Es decir, es la suma de obtener exactamanete 0 defectuosos, más exactamente 1 defectuoso, más exactamnente 2 defectuosos, más exactamente 3 defectuosos.

pbinom(q = 3, size = 10, prob = 0.2)

## [1] 0.8791261

Distribuciones Continuas

Una variable aleatoria X tiene una distribución continua si existe una función f definida en los reales tal que para algún intervalo A
$Pr(X \in A) = \int_{A} f(x)$

La función f es la función de densidad de probabilidad. PDF por sus siglas en inglés.
La probabilidad de que X tome algún valor en un intervalo se encuentra integrando f en ese rango.

Por ejemplo para la distribución uniforme en un intervalo (a,b) podemos ver que su pdf (o función de densidad de probabilidad) es
$f(x) = \begin{cases}\frac{1}{b-a} & \text{para } a \leq x \leq b \\ 0 & \text{resto}\\ \end{cases}$

Distribución Normal

Para la distribución Normal tenemos los siguientes comandos:

Para obtener n variables aleatorias provenientes de una normal con media mean y desvío sd

set.seed(1)
rnorm(n = 5, mean = 10, sd = 2)

## [1]  8.747092 10.367287  8.328743 13.190562 10.659016

Para obtener el valor de la pdf de la normal para algún valor de X en particular. Recuerden que no es una probabilidad, solo es el valor de la función. Las probabilidad se encuentra integrando la función en el intervalo deseado.
Si graficáramos los valores de dnorm para el intervalo -3,3 obtendríamos la forma típica de la normal.

dnorm(0.5, mean = 0, sd = 1)

## [1] 0.3520653

Para obtener la probabilidad acumulada hasta determinado punto. También conocido como Función de Distribución o Función de Distribución Acumulada C.D.F. por sus siglas en ingles Por ejemplo, cual es la probabilidad de obtener un valor igual o menos a 1.5 si tomamos una muestra de una normal estándar

$N \sim (0,1)$

pnorm(q = 1.5, mean = 0, sd = 1)

## [1] 0.9331928

Hay 93.31% de chances de obtener un valor inferior a 1.5 si tomamos una muestra al azar de una normal con media 0 y desvío 1.

La inversa también se puede calcular facilmente en R. Que valor debe tomar la variable aleatoria normal si deseo tenes un 93.31% de chances de obtener un valor menor o igual a ese?

qnorm(p = 0.9331, mean = 0, sd = 1)

## [1] 1.499284

La diferencia respecto al código anterior es porque redondeamos la probabilidad.