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 
# na concatenación (`c`) de dos 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

Declaración de funciones:

# Ejemplo: Fibonacci, que admite dos parámetros:
#  - 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. 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.02582443
median(x)
## [1] -0.01313418
var(x)
## [1] 1.005049
sd(x)
## [1] 1.002521

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.67546 -0.01313 -0.02582  0.66400  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

Una forma de analizar y comparar distribuciones, a veces más acertada que la comparación de funciones de densidad teóricas y empíricas, es mediante Q-Q plots (donde Q viene de cuantil). 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 tienen distribuciones muy parecidas.

2.5.1 Ejemplo 1

Generamos 1000 muestras de una distribución U[0,1] (un pulso cuadrado), y la comparamos con la función de densidad teórica:

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

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

Con esta figura sería arriesgado concluir con certeza que los datos siguen la distribución teórica, dada la variabilidad introducida por la estimación de la función de densidad (parecería que hay datos menores que 0 y mayores de 1). Veamos ahora con 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")

Aquí se puede apreciar que ambos conjuntos tienen un comportamiento similar, y la conclusión de que son la misma variable aleatoria puede emitirse con más certeza.

2.5.2 Ejemplo 2

Generamos ahora dos conjuntos de 1000 valores, de una distribución exponencial y de la suma de tres distribuciones uniformes

x <- rnorm(1000, 0, 1)
y <- runif(1000, -1, 1) + runif(1000, -1, 1) + runif(1000, -1, 1)

Si se comparan sus funciones de distribución, podrían parecer muy similares

plot(ecdf(x))
lines(ecdf(y), col='blue')

Pero el Q-Q plot muestra que en los extremos parece que no hay coincidencia con la recta \(y=x\)

qqplot(x, y)
abline(0, 1, lty=2, col="red")

Para confirmar que estas diferencias no son debidas al azar (en los cuantiles altos y bajos hay menos muestras para poder tener certeza), se puede incrementar la precision si se aumenta el número de muestras:

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

Resulta interesante cambiar la forma en que se generan \(x\) e \(y\) (por ejemplo: aumentando la varianza, o bien sumando una constante) y ver el impacto en el Q-Q plot.

3 Teorema de Bayes (opcional)

La probabilidad condicionada se define como

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

Intercambiando las varibles \(B\) y \(A\) se llega a

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

De donde se obtiene

\[ \Pr(A, 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.01040222

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é?

4 Resolución de problemas con R

4.1 Problema 1.7

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.

1. 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

Para simular este problema en R, primero, 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 ``famosa’’; de hecho, tiene varios modos.

Por el teorema central del límite, la media seguirá una distribución normal: 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 tiene por media un resultado muy parecido al valor teórico (1085)

mean(pkts_avgs)
## [1] 1085.758

Se puede aplicar un test de Student, o t-test, para obtener una estimación de la media real con un intervalo de confianza:

t <- t.test(pkts_avgs)
t
## 
##  One Sample t-test
## 
## data:  pkts_avgs
## t = 2074.3, df = 999, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  1084.731 1086.785
## sample estimates:
## mean of x 
##  1085.758
# valor teórico: 1085
plot(density(pkts_avgs))
abline(v=t$estimate, lty=2, col="red")
abline(v=t$conf.int)

Se deja como ejercicio el siguiente apartado del problema:

2. Calcule la longitud media del tráfico que no es P2P.

4.2 Problema 1.8

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)]\]

La función pmin calcula el mínimo de dos vectores valor a valor. Se puede confirmar que la función de densidad del mínimo de dos uniformes no es sigue una uniforme:

u1 <- runif(1000, 0, 1)
u2 <- runif(1000, 0, 1)

umin <-pmin(u1, u2)
plot(density(umin))

Y parece que tiene por media un valor próximo a 1/3:

mean(umin)
## [1] 0.3502118

Al igual que en el caso anterior, para ser más rigurosos se puede repetir el anterior experimento un número elevado de veces, y así hacer un test de Student para obtener una estimación para la media con su intervalo de confianza:

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 = 1414.5, df = 999, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.3328434 0.3337682
## sample estimates:
## mean of x 
## 0.3333058

Se deja como ejercicio el siguiente problema:

Problema 1.5. Obtenga la expresión de la función de densidad del mínimo de tres variables aleatorias independientes distribuidas uniformemente entre 0 y 1.

4.3 Problema 1.9

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.

Se puede empezar por simular el retardo en la red un número de muestras y estimar su media y función de densidad

N <- 10000
retardo <- sample(c(
    runif(0.40 * N, 10, 70),
    rexp (0.60 * N, 1/30)
  ))
plot(density(retardo))

mean(retardo)
## [1] 33.74754

Pero interesan únicamente las muestras mayores de 50~ms, que es obviamente mayor

retardo50 <- retardo[retardo>50]
plot(density(retardo50))

mean(retardo50)
## [1] 68.9346

Como en casos anteriores, se puede mejorar el rigor de la respuesta a base de repetir el experimento muchas veces

f <- function(i){
  N <- 10000
  retardo <- sample(c(
    runif(0.40 * N, 10, 70),
    rexp (0.60 * N, 1/30)
  ))
  return(mean(retardo[retardo>50]))
}
out <- sapply(1:1000, f)
t <- t.test(out)
t
## 
##  One Sample t-test
## 
## data:  out
## t = 4678.8, df = 999, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  69.16032 69.21835
## sample estimates:
## mean of x 
##  69.18933