2019¶

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-20
in matematica, algebra
7 min read

Esencia del Algebra Lineal

El álgebra lineal está por todas partes en estadística y data science. Matrices, vectores y transformaciones son términos que se escuchan seguido y están detrás de muchos de los métodos y algoritmos que se usan hoy por hoy. Aunque no sea necesario saber del tema para correr un modelo empaquetado en una librería de R, es muy útil entender lo que hacemos realmente ya que nos permite ver a los modelos como algo lógico y no una caja negra mágica.

Hay una serie de videos excelente en inglés que mediante visualizaciones y animaciones permite entender la intuición de muchos de los conceptos básicos, que solo con un libro puede ser medio críptico o poco imaginable. Para el que le interese: ESSENCE OF LINEAR ALGEBRA por 3Blue1Brown.

Este post, aunque quizás medio desordenado y sin mucha prolijidad, es una recopilación de algunas notas. Puede que queden términos en inglés intercalados.

Matrices y vectores

Vector vive en n dimensiones.
Suma de vectores es combinación lineal
En $R^{2}$ $\hat \imath = \left< 1, 0 \right> \text{y} \hat \jmath = \left< 0, 1 \right>$ forman una base. Cualquier punto es una combinación lineal de ellos.
Span es el espacio que pueden generar x vectores. $R^{2}$ es el span de $\left< 1, 0 \right> \left< 0, 1 \right>$
Vector puede ser pensado como una flecha desde el origen (0,0) a las coordenadas que lo identifican. O como un punto directo en las coordenadas..
Matriz es una transformación. Lleva un vector a otro punto. Si transformamos cada posible vector de un espacio por la matriz podemos ver como el espacio es transformado. Ej: rotar, invertir, estirar.
Si transformamos una base, cada punto nuevo puede generarse transformando la nueva base.
Por ej: $z = \left< 3, 2 \right> \text{es } 3\begin{bmatrix} 1 \\ 0 \end{bmatrix} + 2\begin{bmatrix} 0 \\ 1 \end{bmatrix} = 3\hat \imath + 2\hat \jmath$
Aplicando la transformación de la matriz A = $\begin{pmatrix} A & B \\ C & D \end{pmatrix} \text{obtenemos los nuevos vectores base } \hat \imath^{*} \text{y } \hat \jmath^{*}$
$z^{*} = 3\hat \imath^{*} + 2\hat \jmath^{*}$
Multiplicar 2 matrices es transformar un espacio con la primera matriz ( desde la derecha) y luego transformar el resultado por la segunda matriz. Ej: Rotar un espacio y luego invertirlo.
AB != BA -> El orden de las transformaciones importa y se lee de derecha a izquierda.
La matriz (transformación) ya dice como van a ser las nuevas bases.
Si la matriz es $\begin{pmatrix} A & B \\ C & D \end{pmatrix}$, el nuevo $\hat \imath^{*}$ es $\begin{bmatrix} A \\ C \end{bmatrix}$ y $\hat \jmath^{*}$ es $\begin{bmatrix} B \\ D \end{bmatrix}$
Ej: $z = \left< 3, 2 \right> z^{*} = \begin{bmatrix} A & B \\ C & D \end{bmatrix}\begin{bmatrix} 3 \\ 2 \end{bmatrix} = \begin{bmatrix} 3A + 2B \\ 3C + 2D \end{bmatrix}$
Se puede ver también como: $\(z = 3\hat \imath + 2\hat \jmath \text{ } z^{*} = 3\hat \imath^{*} + 2\hat \jmath^{*} = 3\begin{bmatrix} A \\ C \end{bmatrix} + 2\begin{bmatrix} B \\ D \end{bmatrix} = \begin{bmatrix} 3A + 2B \\ 3C + 2D \end{bmatrix}$\)
!!!. Las transformaciones afectan el area (en R2, el volumen en R3..) de las figuras en el espacio (todas por igual). El DETERMINANTE de una matriz mide ese cambio.
Si el determinante es 0 significa que se perdió una dimensión o que todo se comprimió. Pasa de $R^{2}$ a una recta (o a un punto!)
Si el determinante ** es < 0** significa que el espacio se invirtió (en sentido.. como dar vuelta una hoja) pero |DET| siguen siendo el cambio en el area.
A^-1A = I -> una transformación que no hace nada.
Si DET(A) = 0 no existe la matriz inversa. Ej. $R^{2}$ -> si det(A) = 0 la transformación lleva el espacio a una recta. No hay función que lleve cada vector de la recta a un punto en $R{2}$. No hay vuelta atrás.
Si una transformación lleva todos los puntos a una recta tiene rango 1, si lleva a un plano rango 2, y así.. RANGO es el número de dimensiones del output. Rango completo es cuando mantiene las dimensiones del input.
El conjunto de posibles outputs de $A\vec v$ es el Column Space = Span de las columnas
Cuando perdés dimensiones por la transformación todo un conjunto de vectores pasa a ser (0,0). Eso se llama Null Space o Kernel
Matrices no cuadradas cambian la dimensión del espacio.
$$ \begin{bmatrix} A & D \ B & E \ C & F \end{bmatrix} \begin{bmatrix} 1 \ 1 \end{bmatrix} = \begin{bmatrix} A + D \ B + E \ C + F \end{bmatrix} $$ Quedan todos los puntos de $R^{2}$ en un plano en el espacio $R^{3}$. De acá viene la restricción para multiplicar matrices. La cantidad de columnas de la transformación tiene que ser igual a la dimensión del input

DOT PRODUCT o PRODUCTO INTERNO

Dot product entre dos vectores equivale a proyectar uno en el otro y multiplicar sus largos. $\vec A \cdot \vec B = |A^{*}| * |B|$
$A^{*}$ es el vector A proyectado en B.
$\vec B$ es un vector 2D pero también se lo puede ver como una matriz 1x2 que lleva del 2D a la recta.
$\vec B \cdot \vec A = B \vec A \text{que sería llevar A al espacio transformado por B.}$
$B = \begin{bmatrix} B_x & B_y \end{bmatrix}$ tiene en sus columnas donde queda $\hat \imath \text{y } \hat \jmath$ (los vectores unitarios) al ser transformados o algun valor escalado de esto.
$\vec A \cdot \vec B$ es el valor de A en la recta a la que te lleva la transformación B.
Es equivalente proyecto B en A.
Si Dot Product > 0, tienen dirección similar.
Si Dot Product = 0, son ortogonales - proyección que cae en el origen.
Si Dot Product < 0, tienen direcciones opuestas.

CROSS PRODUCT

Está definido para vectores en $R^{3}$
El cross product $\vec u \times \vec v$ es el area del paralelograma que se puede imaginar con las paralelas de los vectores (imaginandolo en $R^{2}$. El signo depende de la orientación de los vectores. El vector de la "derecha" tiene que estar primero para que el cross product sea > 0.
En realidad el paralelograma formado por dos vectores en R³ tiene area equivalente al Largo del vector output de su cross product. Ese vector es ortogonal al paralelograma.

CAMBIO DE BASE

Distintos sistemas de coordenadas definen $\hat \imath = \left< 1, 0 \right>, \hat \jmath \left< 0, 1 \right>$ como algo distinto. NO hay una sola "grilla" válida. El espacio no tiene grilla predeterminada.
Un mismo vector tiene distintas coordenadas según el sistema de bases desde donde se lo mire.
Para pasar de una base a otra se aplica una transformación lineal.
Si $\vec v$ es un vector que queremos pasar de una base a otra, lo transformamos por la nueva base.
Y $\hat \imath^{*} = \left< \hat \imath^{*}_1, \hat \imath^{*}_2 \right>, \hat \jmath^{*} = \left< \hat \jmath^{*}_1, \hat \jmath^{*}_2 \right>$ Entonces: $\(\begin{bmatrix} \hat \imath^{*}_1 & \hat \jmath^{*}_1 \\ \hat \imath^{*}_2 & \hat \jmath^{*}_2 \end{bmatrix} \begin{bmatrix} v_1 \\ v_2 \end{bmatrix} = \begin{bmatrix} v^{*}_1 \\ v^{*}_2 \end{bmatrix}$\)
Donde $\begin{bmatrix} v^{*}_1 \\ v^{*}_2 \end{bmatrix}$ es el vector en la nueva base, es decir, serían las coordenadas del vector $\vec v$ en el nuevo sistema de coordenadas y representando ese punto bajo el sistema de coordenadas original. -> Como se vería $\vec v$ en la nueva base? desde un punto de vista de la base original.
La matriz transforma un vector siguiendo en el lenguaje de la base original.
Ej: Si $\vec v$ es (1,2) en el sistema cartesiano típico y aplicamos la matriz de cambio de base, un vector (1,2) bajo otros ejes se ubicaría en otro punto del espacio. Qué punto es ese bajo el sistema cartesiano? Es (1,2) en el nuevo, pero queremos saber su equivalente en el sistema original.
Por otra parte si queremos saber que coordenadas tomaría el vector $\vec v$ bajo otra base debemos multiplicar por la inversa de la transformación. Transforma el vector al lenguaje de la nueva base. Responde a la pregunta. Qué coordenadas toma el punto V_1, V_2 del espacio en el sistema nuevo?
Para aplicar una transformación a otra base conviene llevar el vector a transformar a la base original, transformar y reconvertir a la nueva base. $$ [A]^{-1}[T][A]\vec v = \vec v^{*}$$
A lo expresa en términos de la base original, luego se le aplica la transformación T y luego se lo devuelve al lenguaje de la nueva transformación.

Eigenvalues y Eigenvectors (autovalores y autovectores)

!!! Al aplicar una transformación lineal a un espacio algunos vectores no cambian de dirección, solo se estiran o contraen pero sobre la misma recta. El resto sí se mueve. Los que se mantienen son los eigenvectors, y su factor de expansión o contracción es su eigenvalue.
Si A es la matriz de transformación, $\vec v$ es un eigenvector y $\lambda$ su eigenvalue. $$ A\vec v = \lambda \vec v$$
Si una transformación es una matriz diagonal, lo único que hace es estirar $\hat \imath \text{y } \hat \jmath$ por lo tanto los vectores base son eigenvectors y la diagonal son los eigenvalues.
Si cambiamos la base a una formada por los eigenvectors (que spanean el espacio) de la matriz podemos aplicar la transformación (la matriz original de donde salieron los eigenvectors) a esta nueva base y solo la va a estirar, por lo tanto es una transformación diagonal, que permite calculos mucho más fácil. Después habría que volver a la base original.
A -> Matriz de transformación E -> Matriz de autovectores que forman la nueva base $\begin{bmatrix} e_11 & e_21 \\ e_12 & e_22 \end{bmatrix}$ D -> Matriz Diagonal cuyos valores son los eigenvalues $$ E^{-1}AE=D$$

E cambia la base a eigenvectors (expresado en la base original), A aplica transformación y E^-1 lo lleva al lenguaje de la nueva base (queda expresado en las nuevas coordenadas)

Espacios Vectoriales Abstractos

!!! Ver funciones como un tipo especial de vectores.
Las funciones se pueden sumar y escalar $f(x) + g(x) \text{y } 2f(x)$
Existen transformaciones lineales de funciones, convierten una función en otra. También conocidas como "operadores"
Para que una transformación sea lineal tiene que cumplir aditividad y mulitplicación por escalar
$$ L(\vec v + \vec w) = L(\vec v) + L(\vec w)$$ $$ L(c\vec v) = cL(\vec v)$$
En general cualquier espacio que cumpla los axiomas los espacios vectoriales puede ser considerado uno.

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.

2019-03-30
in R
6 min read

Introduccion a graficos con mapas

Data

Vamos a ver un ejemplo sencillo para representar información visualmente sobre mapas. En este caso un pequeño dataset de incendios forestales en Argentina de 2012 a 2015. La idea es usar ggplot y mantener el enfoque de gráficos por capas.

Vamos a necesitar.

> tidyverse
> rgdal
> rgeos

Tendremos como input las provincias, departamento, cantidad de focos por incendio, area afectada y año de inicio y fin. Cada observación es un incendio.

Para este ejemplo nos vamos a centrar en las provincias, los focos y su efecto sin importar la fecha.

Empezamos cargando la data.

library(tidyverse)
# Load Raw data
raw <- read.csv("../../static/post/2019-03-30-introduccion-a-graficos-en-mapas/focosincendio.csv", sep = ";")
raw <- as.tibble(raw)

## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## i Please use `as_tibble()` instead.
## i The signature and semantics have changed, see `?as_tibble`.

glimpse(raw)

## Rows: 120
## Columns: 11
## $ pais_id         <int> 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32,~
## $ pais            <chr> "Argentina", "Argentina", "Argentina", "Argentina", "Argentina", "Argentina", "Argentina", "Argentina", "Argentina", "Argentina~
## $ provincia_id    <int> 6, 14, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 30, 30, 30, 30, 30, 30, 30, 30, 30, ~
## $ provincia       <chr> "Buenos Aires", "Córdoba", "Corrientes", "Corrientes", "Corrientes", "Corrientes", "Corrientes", "Corrientes", "Corrientes", "C~
## $ departamento_id <int> 833, 14, 56, 28, 28, 28, 70, 84, 84, 84, 112, 112, 119, 119, 126, 147, 154, 154, 161, 168, 168, 168, 168, 8, 15, 15, 15, 15, 28~
## $ departamento    <chr> "Tres Arroyos", "Calamuchita", "General Alvear", "Concepción", "Concepción", "Concepción", "Goya", "Itizaingó", "Ituzaingo", "I~
## $ sup_afectada    <dbl> 2400.00, 50.00, 257.00, 130.00, 5.00, 146.00, 30.00, 294.30, 378.00, 158.00, 300.00, 450.00, 450.00, 15.00, 20.00, 141.00, 295.~
## $ uni_med_id      <chr> "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "~
## $ cant_focos      <int> 1, 1, 1, 2, 1, 1, 1, 3, 1, 3, 1, 1, 5, 2, 2, 1, 3, 1, 1, 18, 2, 3, 2, 1, 0, 0, 1, 7, 0, 2, 1, 0, 0, 0, 5, 10, 10, 3, 1, 1, 1, 1~
## $ año_inicial     <int> 2014, 2015, 2012, 2012, 2013, 2015, 2012, 2012, 2013, 2014, 2012, 2013, 2012, 2013, 2015, 2012, 2012, 2014, 2012, 2012, 2013, 2~
## $ año_final       <int> 2014, 2015, 2012, 2012, 2013, 2015, 2012, 2012, 2013, 2014, 2012, 2013, 2012, 2013, 2015, 2012, 2012, 2014, 2012, 2012, 2013, 2~

Exploramos un poco el dataset. Lo que nos vas a interesar represetnar es la segunda parte del código. Variables agregadas a nivel provincia.

# Generate Summary to Explore
sum_year <- raw %>% group_by(año_inicial) %>%
  summarise(focos = sum(cant_focos), sup = sum(sup_afectada, na.rm = TRUE)) %>%
  mutate(sup_prom = sup/focos)

# Actual data to be plotted
sum_prov <- raw %>% group_by(provincia) %>%
  summarise(focos = sum(cant_focos), sup = sum(sup_afectada, na.rm = TRUE)) %>%
  mutate(sup_prom = sup/focos) %>%
  arrange(desc(sup_prom)) %>%
  mutate(provincia = as.character(provincia))

head(sum_prov)

## # A tibble: 6 x 4
##   provincia    focos   sup sup_prom
##   <chr>        <int> <dbl>    <dbl>
## 1 Buenos Aires     1 2400    2400  
## 2 Corrientes      55 5440.     98.9
## 3 Misiones        10  679      67.9
## 4 Córdoba          1   50      50  
## 5 Entre Ríos      12  594      49.5
## 6 Jujuy           62 2721.     43.9

Tenemos por provincia la cantidad de focos y la superficie afectada. Vamos a usar la superficie promedio por foco para visualizar la magnitud de cada foco de incendio.

Ahora llega lo importante. ¿Cómo representar esta data provincial en un mapa de Argentina?

Mapa

Necesitamos un ShapeFile de Argentina, que basicamente es el tipo de archivo que se usa para representar mapas en gráficos. Contiene divisiones del país (provincias) con sus respectivas coordenadas y nombres.

Utilizaremos data descargada del siguiente link. url <- "http://biogeo.ucdavis.edu/data/diva/adm/ARG_adm.zip"

Lo descargué y deszipee en la computadora. Lo leemos con una librería particular RGDAL.

dsn contiene la ruta a la carpeta con los archivos del shapefile. Layer apunta a al set de archivos que contiene la data que queremos. Generalmente hay otros sets con información no relevante al gráfico.

argentina <- rgdal::readOGR(dsn = "../../static/post/2019-03-30-introduccion-a-graficos-en-mapas/ARG_adm", layer = "ARG_adm1", use_iconv=TRUE, encoding='UTF-8')

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\DataScience\StatsBlog\blogStats\static\post\2019-03-30-introduccion-a-graficos-en-mapas\ARG_adm", layer: "ARG_adm1"
## with 24 features
## It has 9 fields
## Integer64 fields read as strings:  ID_0 ID_1

Es un archivo S4 por lo que se utiliza "@" para acceder a su contenido. Por ejemplo:

head(argentina@data)

##   ID_0 ISO    NAME_0 ID_1                 NAME_1           TYPE_1        ENGTYPE_1 NL_NAME_1
## 0   12 ARG Argentina    1           Buenos Aires        Provincia         Province      <NA>
## 1   12 ARG Argentina    2                Córdoba        Provincia         Province      <NA>
## 2   12 ARG Argentina    3              Catamarca        Provincia         Province      <NA>
## 3   12 ARG Argentina    4                  Chaco        Provincia         Province      <NA>
## 4   12 ARG Argentina    5                 Chubut        Provincia         Province      <NA>
## 5   12 ARG Argentina    6 Ciudad de Buenos Aires Distrito Federal Federal District      <NA>
##                                                                               VARNAME_1
## 0                                                                   Baires|Buenos Ayres
## 1                                                                               Cordova
## 2                                                                                  <NA>
## 3                                                        El Chaco|Presidente Juan Peron
## 4                                                                                  <NA>
## 5 BUENOS AIRES D.F.|Capital Federal|Distretto Federale|Distrito Federal|Federal Capital

Este tipo de archivos tiene una estructura complicada y hay varias librerías útiles. Con el fin de mantenernos dentro del tidyverse usaremos el enfoque de ggplot por capaz para graficar. Primero necesitamos llevar la información del shapefile a un dataframe.

# Transformo a dataframe.
argentina_df <- broom::tidy(argentina)

## Regions defined for each Polygons

# id es la provincia
head(argentina_df)

## # A tibble: 6 x 7
##    long   lat order hole  piece group id   
##   <dbl> <dbl> <int> <lgl> <fct> <fct> <chr>
## 1 -60.2 -33.3     1 FALSE 1     0.1   0    
## 2 -60.2 -33.3     2 FALSE 1     0.1   0    
## 3 -60.2 -33.3     3 FALSE 1     0.1   0    
## 4 -60.2 -33.3     4 FALSE 1     0.1   0    
## 5 -60.2 -33.3     5 FALSE 1     0.1   0    
## 6 -60.2 -33.3     6 FALSE 1     0.1   0

Ahora genero un diccionario de ids con su respectiva provincia para poder linkear mi mapa con la data de incendios.

ids <- cbind.data.frame(provincia = as.character(argentina@data$NAME_1), id = as.character(rownames(argentina@data))) %>%
  mutate(provincia = as.character(provincia), id = as.character(id))

El dataframe generado a partir del shapefile tiene las coordenadas que forman cada provincia y el id. Ahora lo que haremos es pegarle la data de incendios para poder utilizarla sobre el mapa. Está generado el log de los focos porque en una prueba intenté usar la data transformada para suavizar outliers pero no se verá en esta versión.

# Agrego provincia por ID y data de incendios a la data del shapefile.
argentina_df2 <- argentina_df %>% left_join(ids, by = "id") %>%
  left_join(sum_prov, by = "provincia") %>%
  mutate(focos = ifelse(is.na(focos),0.5,focos),
         logfocos = log(focos), 
         logfocos2 = logfocos - min(logfocos))

Por ahora vamos a poder ser capaces de graficar un mapa de argentina con sus provincias delimitadas y pintarlas según la cantidad de focos de incendio por ejemplo. Pero además vamos a querer agregar alguna forma sobre cada provincia. Por ejemplo un punto de distinto tamaño según el area afectada por cada foco en promedio. Para eso necesitamos localizar el centroide de cada poligono, es decir, el centro de cada provincia. Para ellos usamos la librería RGEOS Luego a cada centro le agrego la data que voy a querer usar. En este caso la superficie promedio afectada por foco de cada provincia.

# Calculo el centro de cada poligono (provincias)
# para obtener el "centro" donde iran los puntos o nombres.
centros <- rgeos::gCentroid(argentina, byid = TRUE) %>%
  as.data.frame() %>%
  mutate(id = rownames(.))

# Agrego data relevante para el ploteo (superficie promedio)
centros2 <- centros %>% left_join(ids, by = "id") %>%
  left_join(sum_prov, by = "provincia") %>%
  mutate(focos = ifelse(is.na(focos),0.5,focos),
         sup = ifelse(is.na(sup),0,sup),
         log_sup_prom = log(sup_prom),
         sup_prom_sin_outlier = ifelse(sup_prom > 150, 150,sup_prom) ) # esto es para suavizar el outlier. Buscar otro enfoque.

Grafico!

Ya tenemos todo. Tenemos el mapa, tenemos la cantidad de focos de incendio por provincia y tenemos el centro de cada provincia donde vamos a incluir un punto que muestra la intensidad de los incendios. Simplemente graficamos sigueindo la lógica por capas de ggplot. Las provincias sin puntos son aquellas que no tuvieron ningún incendio.

ggplot() +
  geom_polygon(data = argentina_df2, aes(x=long, y = lat, group = group,  fill = focos), color = "white") + # mapa de argentina
  # coloreado segun cantidad de focos
  coord_fixed(0.8) + # tamañp del mapa
  scale_fill_gradient2("Cantidad de Focos de incendio", low = "white", mid = "lightgreen", high = "darkred") + # escala de colores para focos
  geom_point(data = centros2, aes(x = x, y = y, size = sup_prom_sin_outlier)) + # puntos por provincia con superficie promedio
  scale_size(name = "Superficie Promedio Afectada (ha)",range = c(1,5)) + # escala de los puntos
  guides(fill = guide_legend(order = 1), # Orden de los leyendas a la derecha.
         size = guide_legend(order = 2)) # Por algun motivo esto discretizo la leyenda de focos