Licenciado bajo Creative Commons BY-NC-SA 4.0.

1 Introducción a R

1.1 ¿Qué es R?

R es un lenguaje de programación de scripting orientado al análisis estadístico y gráfico. Cuenta con una nutrida masa de usuarios y desarrolladores que lo han convertido en líder indiscutible dentro de la comunidad estadística internacional, tanto en el ámbito empresarial como en el académico.

1.2 Entorno de programación

Vamos a utilizar R desde la aplicación gráfica RStudio, un IDE (Integrated Development Environment) que proporciona un espacio de trabajo muy amigable y muy similar al de MATLAB. Al abrir el programa, se distinguen cuatro regiones de arriba a abajo y de izquierda a derecha:

Editor de scripts. Se recomienda trabajar aquí para guardar un archivo coherente con los ejercicios que se vayan realizando. Al seleccionar código y darle a Run, automáticamente se ejecuta en el panel de debajo, la consola.
Consola de comandos, para introducción directa de código o ejecución desde el editor.
Panel de historial de comandos y workspace similar al de MATLAB, donde se visualizan las variables en uso y su valor (escalares) o tipo (vectores, matrices, etc.).
Panel en el que se agrupan las pestañas de gestión de archivos, gráficas, paquetes y ayuda.

1.3 R basics

# Esto es un comentario

# Un cálculo simple
2^10 - 24

## [1] 1000

Operador de asignación: <- (= también es válido)

x <- 1
x = 1

Si se quiere poner dos instrucciones en la misma línea, hay que separar con punto y coma. Escribiendo solo la variable en la terminal mostraremos su valor por terminal

# Two commands in the same line
x <- 2^10 - 24; x

## [1] 1000

Creación de secuencias:

# secuencia simple
x <- seq(1, 2, 0.1); x

##  [1] 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0

# secuencias más compleja usando la función seq(inicio, final, paso)
y <- seq(1, 2, 0.1); y

##  [1] 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0

#Una secuencia algo más heterogenea, hecha como una concatenación (`c`) de partes (1, 5 y seq(100, 110, 2))
z <- c(1, 5, seq(100, 110, 2)); z

## [1]   1   5 100 102 104 106 108 110

Acceso a elementos de una secuencia:

x <- 1:10;

# Accedemos al primer elemento de la secuencia
x[1]

## [1] 1

¿Qué hace tal o cuál función? La ayuda de R es muy completa. Tenemos la búsqueda específica por nombre de función:

? c

## starting httpd help server ... done

Declaración de funciones:

#  - n:     numero de valores de la secuencia que queremos calcular (required)
#  - full:  Indica si quermos que la función devuelva toda la secuencia o solo el último valor (optional)  DEFAULT=FALSE
fibonacci <- function(n, full=FALSE) {
  # Reservamos memoria para el vector
  fibvals <- numeric(n)
  
  # ponemos los dos primeros valores de la secuencia
  fibvals[1] <- 1
  fibvals[2] <- 1
  
  # Bucle para calcular el resto de valores
  if (n > 2) {
    x <- 3:n
    for (i in x) { 
      fibvals[i] <- fibvals[i-1] + fibvals[i-2]
    }
    
  }
  
  # Por defecto, devolvemos el valor en la posicion n
  if (!full) {
    return(fibvals[n])
  } 
  # en otro caso, retornamos toda la secuencia
  else {
    return(fibvals)
  }
}

# Pedimos el décimo elemento
fibonacci(10)

## [1] 55

 # Toda la serie hasta el décimo elemento
fibonacci(10, TRUE)

##  [1]  1  1  2  3  5  8 13 21 34 55

2 Estadística con R

2.1 Generación de variables aleatorias

Los paquetes base y stats son fundamentales en R (vienen por defecto) y nos proveen de gran cantidad de funciones básicas que en otros lenguajes deberíamos implementar por nuestra cuenta. Echemos un vistazo a las distribuciones disponibles. Léase lo siguiente, donde disponemos de links a la ayuda de las diferentes distribuciones:

? distributions

Vemos que, para cada distribución xxx, hay cuatro funciones disponibles:

rxxx(n, ...) genera \(n\) muestras aleatorias distribuidas según xxx.
dxxx(x, ...) devuelve \(f(x)\) para un valor o conjunto de valores, donde \(f\) es la función de densidad que hemos visto en teor??a (\(Pr\) para distribuciones discretas). Por tanto, para cualquier xxx, dxxx(Inf) == dxxx(-Inf) == 0.
pxxx(q, ...) devuelve \(F(q) = \int_{-\infty}^q f(x)dx\), donde \(F\) es la función de distribución que hemos visto en teoría (basta cambiar la notación integral por sumatorios en el caso de distribuciones discretas). Por tanto, pxxx(Inf) == 1.
qxxx(p, ...) devuelve el \(q\) para el que \(F(q) = p\). Es decir, la función inversa a la anterior, y pxxx(1) == Inf.

Generamos 1000 muestras de una \(N(0,1)\) (normal de media 0, varianza 1) y simplemente vemos lo generado

# Generamos valores random de una distribucion normal
x <- rnorm(1000, 0, 1)

# Mostramos los valores obtenidos
plot(x)

2.2 Estadística descriptiva

Comprobamos media, mediana, varianza y desviación estándar:

mean(x)

## [1] -0.02482999

median(x)

## [1] -0.01313418

var(x)

## [1] 1.009967

sd(x)

## [1] 1.004971

La función summaryes muy versátil, ya que acepta muchos tipos de datos de entrada y da mucha información:

summary(x)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -3.37174 -0.68834 -0.01313 -0.02483  0.66443  3.49530

2.3 Función de densidad

Representa como una variable aleatoria se distribuye. Por ejemplo una distribución \(N(0,1)\) tiene una función de densidad \(f(x)\) con forma de campana de Gauss centrada en 0. Vemos la función de densidad empírica que tienen nuestros datos:

# Empírica
## Generamos datos desde una distribución normal
x <- rnorm(1000, 0, 1)
## Representamos la funcion de densidad de los datos generados
plot(density(x))

# Teórica
## Ahora superponemos la función de densidad teórica de una distribución N(0,1) usando dnorm
##    - add=TRUE  superpone la línea en la ultima imagen generada (en lugar de generar una nueva
##    - col="red" color de la línea
##    - lty=2     tipo de línea
curve(dnorm(x, 0, 1), add=TRUE, lty=2, col="red")

2.4 Función de distribución

De la misma forma, podemos ver la función de distribución F(x) (también llamada CDF) empírica de nuestros datos:

# Empírica CDF
x <- rnorm(1000, 0, 1)
Fx <- ecdf(x)
plot(x, Fx(x))

# Superponemos la teórica usando la función pnorm()
curve(pnorm(x, 0, 1), add=TRUE, lty=2, col="red")

2.5 Q-Q Plots

Otra forma muy utilizada de analizar y comparar distribuciones es mediante la llamada Q-Q plot (donde Q viene de cuantil, y recordemos que estos son el dominio de la función de densidad). Mediante la función qqplot podemos enfrentar dos distribuciones de datos. Si el resultado es aproximadamente una recta de offset 0 y pendiente 1, significa que ambos conjuntos de muestras provienen de la misma distribución

2.5.1 Ejemplo 1

Generamos 1000 muestras de una distribución N(0,1) (pulso cuadrado) y comparamos las funciones de distribución y densidad:

# Empírica
x <- runif(1000, 0, 1)
plot(density(x))

# Teórica
curve(dunif(x, 0, 1), add=TRUE, lty=2, col="red")

Mirando la figura, es dificil concluir que nuestros datos siguen a la distribución teórica. Veamos ahora el Q-Q Plot:

# Generamos el Q-Q plot
qqplot(x, runif(length(x), 0, 1))

# Superponemos una linea con cero offset y pendiente 1
abline(0, 1, lty=2, col="red")

Aqui podemos apreciar que ambos crecen de una forma sismilar. Por esto, podemos decir que ambos conjuntos de datos provienen de la misma variable aleatoria.

2.5.2 Ejemplo 2

Generamos ahora dos conjuntos de 1000 valores de una distribución exponencial y los comparamos con Q-Q plot.

x <- rexp(1000)
y <- rexp(1000)
qqplot(x, y)
abline(0, 1, lty=2, col="red")

Cuando comparamos con Q-Q plot, algunas veces podemos ver “ruido” alrededor de los vertices debido a la falta de muestras en esas regiones (siendo por tanto los quantiles menos precisos). Para incrementar la precision, necesitaremos aumentar el número de muestras:

# Mayor número de muestras
x <- rexp(1000000)
y <- rexp(1000000)
qqplot(x, y)
abline(0, 1, lty=2, col="red")

A la hora de explorar un conjunto de datos, ninguna de ellas es suficiente por separado. Para ilustrar esto, tomemos el conjunto sintético de datos Datasaurus Dozen, disponible en el CRAN:

if (!requireNamespace("datasauRus", quietly=TRUE))
  install.packages("datasauRus")
library(datasauRus)

unique(datasaurus_dozen$dataset)

##  [1] "dino"       "away"       "h_lines"    "v_lines"    "x_shape"   
##  [6] "star"       "high_lines" "dots"       "circle"     "bullseye"  
## [11] "slant_up"   "slant_down" "wide_lines"

head(datasaurus_dozen)

##   dataset       x       y
## 1    dino 55.3846 97.1795
## 2    dino 51.5385 96.0256
## 3    dino 46.1538 94.4872
## 4    dino 42.8205 91.4103
## 5    dino 40.7692 88.3333
## 6    dino 38.7179 84.8718

Tenemos 13 datasets (identificados por la columna dataset), y cada dataset se compone por valores de \(x\) e \(y\). Vamos a calcular las medias y desviaciones estándar de ambas coordenadas para cada dataset:

for (name in unique(datasaurus_dozen$dataset)) {
  dataset <- subset(datasaurus_dozen, dataset==name)
  print(sprintf(
    "mean_x = %.3f, sd_x = %.3f, mean_y = %.3f, sd_y = %.3f  <- %s", 
    mean(dataset$x), sd(dataset$x), mean(dataset$y), sd(dataset$y), name))
}

## [1] "mean_x = 54.263, sd_x = 16.765, mean_y = 47.832, sd_y = 26.935  <- dino"
## [1] "mean_x = 54.266, sd_x = 16.770, mean_y = 47.835, sd_y = 26.940  <- away"
## [1] "mean_x = 54.261, sd_x = 16.766, mean_y = 47.830, sd_y = 26.940  <- h_lines"
## [1] "mean_x = 54.270, sd_x = 16.770, mean_y = 47.837, sd_y = 26.938  <- v_lines"
## [1] "mean_x = 54.260, sd_x = 16.770, mean_y = 47.840, sd_y = 26.930  <- x_shape"
## [1] "mean_x = 54.267, sd_x = 16.769, mean_y = 47.840, sd_y = 26.930  <- star"
## [1] "mean_x = 54.269, sd_x = 16.767, mean_y = 47.835, sd_y = 26.940  <- high_lines"
## [1] "mean_x = 54.260, sd_x = 16.768, mean_y = 47.840, sd_y = 26.930  <- dots"
## [1] "mean_x = 54.267, sd_x = 16.760, mean_y = 47.838, sd_y = 26.930  <- circle"
## [1] "mean_x = 54.269, sd_x = 16.769, mean_y = 47.831, sd_y = 26.936  <- bullseye"
## [1] "mean_x = 54.266, sd_x = 16.769, mean_y = 47.831, sd_y = 26.939  <- slant_up"
## [1] "mean_x = 54.268, sd_x = 16.767, mean_y = 47.836, sd_y = 26.936  <- slant_down"
## [1] "mean_x = 54.267, sd_x = 16.770, mean_y = 47.832, sd_y = 26.938  <- wide_lines"

que, como se ve, son prácticamente iguales hasta 3 cifras decimales. Pero, ¿se parecen entre sí? las distribuciones? En esta ocasión, basta realizar una simple gráfica de cada dataset para verificar que son muy diferentes:

library(ggplot2)

ggplot(datasaurus_dozen, aes(x=x, y=y)) +
  geom_point()+
  theme_minimal()+
  facet_wrap(~dataset, ncol=4)

3 Suma de variables aleatorias independientes

Distribución de la suma de variables aleatorias independientes

# Empírica
x <- rnorm(1000, 0, 1)
y <- rnorm(1000, 10, 3)

z <- x + y
# Mostramos la función de densidad de la suma
plot(density(z))


# Para confirmar que la suma resulta en una convolución, vemos a comparar lo anterior con el resultado de la suma de forma teórica

# Funcion de densidad de x
fx <- function(x){
  dnorm(x, 0, 1)
}
# Función de densidad de y
fy <- function(x){
  dnorm(x, 10, 3)
}
# Función de densidad de la suma
fadd <- function(x, y){
  return(fx(y-x)*fy(x))
}
fz <- function(x){
  # Coge un valor numérico y devuelve un vector numérico del mismo tamaño
  value <- integrate(fadd, -Inf, Inf, x)$value
  return(value)
}
# Vectorizamos la funcion que calcula la integral (para que nuestra funcion esté alineadad con la salida de la integral - vector -
fz <- Vectorize(fz)

curve(fz(x), add=TRUE, lty=2, col="red")

Ahora lo repetimos generamos dos distribuciones uniformes entre 0 y 1, es decir, dos pulsos cuadrados y los sumamos. Sabemos que el resultado de convolucionar dos pulsos cuadrados debe ser un pulso triangular de anchura doble:

# Empírica
x <- runif(1000, 0, 1)
y <- runif(1000, 0, 1)
z <- x + y
plot(density(z))

# Teórica
fx <- function(x){ 
  dunif(x, 0, 1)
}
fy <- function(x){ 
  dunif(x, 0, 1)
}
fadd <- function(x, y){
  return(fx(y-x)*fy(x))
}
fz <- function(x){ 
  value <- integrate(fadd, -Inf, Inf, x)$value
  return(value)
}
fz <- Vectorize(fz)
curve(fz(x), add=TRUE, lty=2, col="red")

4 Teorema de Bayes

En clases de teoría hemos visto cómo calcular la probabilidad condicionada:

\[\Pr(A | B) = \frac{\Pr(A \cap B)}{\Pr(B)}\]

Pero la fórmula no distingue entre variables \(A\) o \(B\), por lo que podríamos escribir igualmente

\[\Pr(B | A) = \frac{\Pr(A \cap B)}{\Pr(A)}\]

Mediante las fórmulas anteriores, obtenemos entonces que

\[ \Pr(A \cap B) = \Pr(A | B)\Pr(B) = \Pr(B | A)\Pr(A)\]

Y de aquí se extrae el famoso Teorema de Bayes:

\[ \Pr(A | B) = \frac{\Pr(B | A)\Pr(A)}{\Pr(B)}\]

Bayes’s theorem is to the theory of probability what Pythagoras’s theorem is to geometry. (Sir Harold Jeffreys)

Vamos a ilustrar la utilidad de tan valioso teorema con un pequeño ejemplo. Supongamos que se ha desarrollado un test para detectar cierta enfermedad incurable cuya prevalencia es del 1/10000, es decir, 1 de cada 10000 personas desarrolla dicha enfermedad. El test tiene una sensibilidad (probabilidad de detección) del 99 % (es decir, de cada 100 personas enfermas, da positivo en 99) y una especificidad del 99 % (de cada 100 personas sanas, da negativo en 99).

Ahora bien, nos hacen el test y da positivo. ¿Qué probabilidad tenemos de tener la enfermedad? Pongamos las probabilidades que conocemos hasta el momento en la notación apropiada:

Conocemos la probabilidad de desarrollar la enfermedad: \(\Pr(\mathrm{enfermo}) = 1/10000 = 0.0001\).
La sensibilidad es la probabilidad de que el test sea positivo condicionada a que el paciente esté enfermo: \(\Pr(+ | \mathrm{enfermo}) = 0.99\).
La especificidad es la probabilidad de que el test sea negativo condicionada a que el paciente esté sano: \(\Pr(- | \mathrm{sano}) = 0.99\).

El objetivo del problema es calcular la probabilidad que tenemos de tener la enfermedad condicionada a que el test ha dado positivo. El Teorema de Bayes nos dice lo siguiente:

\[ \Pr(\mathrm{enfermo} | +) = \frac{\Pr(+ | \mathrm{enfermo})\Pr(\mathrm{enfermo})}{\Pr(+)}\]

Todavía no conocemos la probabilidad de que el test salga positivo, pero podemos calcularla fácilmente sumando los positivos dados por la sensibilidad y los dados por la especificidad:

\[\begin{aligned} \Pr(+) &= \Pr(+ | \mathrm{enfermo})\Pr(\mathrm{enfermo}) + \Pr(+ | \mathrm{sano})\Pr(\mathrm{sano}) \\ &= \Pr(+ | \mathrm{enfermo})\Pr(\mathrm{enfermo}) + \left(1-\Pr(- | \mathrm{sano})\right)\left(1-\Pr(\mathrm{enfermo})\right) \\ & = 0.99\cdot0.0001 + (1-0.99)\cdot(1-0.0001) = 0.010098 \end{aligned}\]

Por último, el resultado buscado:

\[ \Pr(\mathrm{enfermo} | +) = \frac{0.99\cdot0.0001}{0.010098} \approx 0.01 \]

¿Tan solo un 1 % de probabilidad de tener la enfermedad a pesar de haber dado positivo en un test tan aparentemente bueno? Bayes no miente, pero nosotros somos unos escépticos, así que lo vamos a simular con la ayuda de R.

# Funcion para generar datos de personas enfermas o sanas
population <- function(n) {
  # Vector (TRUE/FALSE)
  sick <- runif(n) < 0.0001
  people_df <- data.frame(sick)
  return(people_df)
}

# Ahora, creamos una funcion que simula el resultado de un test con  los datos descritos anteriormente
test <- function(people_df, sensitivity, specificity) {
  random <- runif(nrow(people_df))
  people_df$testedPositive <- (people_df$sick & (random < sensitivity)) |
                           (!people_df$sick & (random > specificity))
  return(people_df)
}

# Generamos una población con tamaño similar a la población de Madrid
madrid <- population(3000000)
# El resultado es un dataframe de una columna que almacena quien está enfermo
head(madrid)

##    sick
## 1 FALSE
## 2 FALSE
## 3 FALSE
## 4 FALSE
## 5 FALSE
## 6 FALSE

# Todos han pasado el test
madrid <- test(madrid, 0.99, 0.99)
# Indicamos en una nueva columna si el test resultó positivo para el individuo
head(madrid)

##    sick testedPositive
## 1 FALSE          FALSE
## 2 FALSE          FALSE
## 3 FALSE          FALSE
## 4 FALSE          FALSE
## 5 FALSE          FALSE
## 6 FALSE          FALSE

# Ahora, calculamos la cantidad de individuos con la enfermedad que dieron positivo en el test
positive <- subset(madrid, testedPositive)
positiveAndSick <- subset(madrid, testedPositive & sick)
nrow(positiveAndSick) / nrow(positive)

## [1] 0.009713275

Y este es el motivo por el que estos tests se repiten.

Queda como ejercicio repetir el test múltiples veces para aquellos que dan positivo. ¿Qué porcentaje de enfermos sobre positivos se obtienen tras un segundo test positivo? ¿Y tras un tercero? ¿Qué es más importante de cara a fiarnos más o menos del resultado final del test, la sensibilidad o la especificidad? ¿Por qué?

5 Resolución de problemas con R

5.1 Problema 1

Suponga que el tráfico en una red se distribuye según se muestra en la tabla, donde la primera columna indica la aplicación, la segunda indica la distribución que siguen sus tramas, y la tercera, la proporción relativa de dichas tramas. Obtenga la longitud media de las tramas en la red.

Application	Length (B)	%
Skype	\(U(50,150)\)	5
P2P	\(U(1000, 1500)\)	60
Web	\(exp(1/1000)\)	25
email	\(N(800,100)\)	10

Adicionalmente, nos podríamos preguntar, ¿se parece el proceso resultante a alguna de las distribuciones anteriores?

Es muy sencillo simular este problema en R. Para ello, generamos un número elevado de muestras de cada proceso, en la proporción adecuada, y las mezclamos aleatoriamente con sample:

# Número total de muestras a generar
N <- 1000
# Generamos la proporción adecuada de cada tipo de tráfico
pkts <- sample(c(
  runif(0.05 * N, 50, 150),
  runif(0.60 * N, 1000, 1500),
  rexp (0.25 * N, 1/1000),
  rnorm(0.10 * N, 800, 100)
))

El resultado es el siguiente, con la media resaltada en rojo:

plot(density(pkts))
abline(v=mean(pkts), lty=2, col="red")

Esta distribución no se parece a ninguna conocida; de hecho, tiene varios modos. Sin embargo, la teoría nos dice que la media empírica seguirá una distribución normal. En otras palabras, si repetimos el proceso muchas veces (generar un nuevo vector de muestras y hallar su media), las medias que vayamos obteniendo seguirán una distribución gaussiana en torno a la media real que queremos obtener.

Esto se podría hacer con un bucle for, pero R proporciona herramientas para sustituir bucles por funciones que trabajan con vectores. En este caso, utilizaremos sapply, que, como su documentación indica, aplica una función sobre una lista o vector, y devuelve un vector de resultados. As?? que introducimos lo anterior en una función para aplicarla 1000 veces:

# Calcula la media de la muestra
gen <- function(i) {
  pkts <- sample(c(
    runif(0.05 * N, 50, 150),
    runif(0.60 * N, 1000, 1500),
    rexp (0.25 * N, 1/1000),
    rnorm(0.10 * N, 800, 100)
  ))
  return(mean(pkts))
}
# El método sapply llama 1000 veces a la función gen() almacenando los resultados en un vector
pkts_avgs <- sapply(1:1000, gen)

Esta variable pkts_avgs contiene las medias de 1000 realizaciones de la simulación inicial. Podemos comprobar que efectivamente estas medias se distribuyen según una gaussiana:

plot(density(pkts_avgs))

Y dado que se cumple esta propiedad, podemos aplicar un test de Student, o t-test, para obtener una estimación de la media real con un intervalo de confianza:

5.2 Problema 2

Sean dos variables aleatorias independientes \(u1\) y \(u2\), uniformemente distribuidas entre 0 y 1. Calcule la esperanza del mínimo de las dos.

\[E[min(\mu_1, \mu_2)]\]

f <- function(i){
  u1 <- runif(1000, 0, 1)
  u2 <- runif(1000, 0, 1)
  return(mean(pmin(u1, u2)))  
}
out <- sapply(1:1000, f)
t <- t.test(out)
t

## 
##  One Sample t-test
## 
## data:  out
## t = 1435, df = 999, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.3325567 0.3334675
## sample estimates:
## mean of x 
## 0.3330121

5.3 Problema 3

El 40% de los paquetes sufren un retardo de red que se puede modelar según una variable aleatoria uniformemente distribuida entre 10~ms y 70~ms, y el 60% restante sufren un retardo según una variable aleatoria exponencial de media 30~ms. Calcule el retardo medio de los paquetes con más de 50~ms de retardo.

delay_threshold <- 50
lambda <- 1/30
get_avg <- function(x){
  u1 <- runif(4000, 10, 70)
  e1 <- rexp(6000, lambda)
  u1 <- u1[u1>50]
  e1 <- e1[e1>50]
  # Combinamos y aleatorizamos el vector
  samples <- sample(c(u1, e1))
  return(mean(samples))
}

average <- sapply(1:1000, get_avg)
mean(average)

## [1] 69.20544

Teoría de Redes - Práctica 1

Iñaki Úcar, Pablo Serrano, José Alberto Hernández, Ginés García Avilés

24 de Enero de 2020