Licenciado bajo Creative Commons BY-NC-SA 4.0.

1 Variable aleatoria exponencial

1.1 Recordatorio de estadística descriptiva

Recordemos brevemente las herramientas básicas con las que trabajamos aplicadas a la variable exponencial. Generamos 1000 muestras de una \(\mathrm{Exp}(5)\) (exponencial de media 5). Tenga en cuenta que rexp acepta como parámetro el “rate”, es decir, “lambda” (la inversa de la media):

x <- rexp(1000, rate=1/2) # rate = lambda = 1/media
plot(x)

summary(x)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##  0.000562  0.668708  1.510648  2.155886  2.917028 16.924672

Función de densidad \(f(x)\) empírica:

plot(density(x), ylim=c(0,0.5))

# f(x) teórica
curve(dexp(x, 1/2), add=TRUE, lty=2, col="red")

Función de distribución \(F(x)\), o CDF, empírica:

# Fx es una función que representa la función de distribución empírica
Fx <- ecdf(x)
# Devuelve la función de distribución empícia en el punto cero (Debería ser 0 porque nuestras muestras vienen de una variable aleatoria exponencial)
Fx(0)
## [1] 0
plot(x, Fx(x))

# F(x) teórica
curve(pexp(x, 1/2), add=TRUE, lty=2, col="red")

En el caso de la variable exponencial, resulta muy útil la función de supervivencia o CDF complementaria: la CCDF (Complementary CDF) o \(F^C(x)\) en la notación usada en clases de teoría. Esto es así porque

\[F^C(x) = 1 - F(x) = 1 - \left[1-e^{-\lambda x}\right] = e^{-\lambda x}\]

Se realiza un logaritmo en la igualdad anterior y se obtiene una recta cuya pendiente es la tasa de la exponencial:

\[\ln F^C(x) = -\lambda x\]

Esto es, mediante comandos con R se puede representar la función de supervivencia o CCDF:

plot(x, 1-Fx(x))

que se convierte en una recta al aplicar el logaritmo:

plot(x, log(1-Fx(x)))

Y de aquí podemos se puede obtener la tasa \(\lambda\) por regresión lineal. (¿Qué es una regresión lineal? Aquí una descripción y un vídeo instructivo).

logCCDF <- log(1 - Fx(x))
# Nótese que la CDF empírica le asigna un 1 a la última muestra de la
# gráfica debido al carácter acumulativo. Eso quiere decir que la CCDF
# tiene un 0 y, al aplicar el logaritmo, obtenemos un -infinito que
# impide  hacer la regresión y que se ha de suprimir tanto en el
# vector de las y como en el de las x
x <- x[is.finite(logCCDF)]
logCCDF <- logCCDF[is.finite(logCCDF)]

# Regresión lineal y resultado
fit <- lm(logCCDF ~ x)
summary(fit)
## 
## Call:
## lm(formula = logCCDF ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.49409 -0.01349  0.00589  0.01265  0.60599 
## 
## Coefficients:
##               Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)  0.0073921  0.0016227    4.555 5.87e-06 ***
## x           -0.4689248  0.0005412 -866.486  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0359 on 997 degrees of freedom
## Multiple R-squared:  0.9987, Adjusted R-squared:  0.9987 
## F-statistic: 7.508e+05 on 1 and 997 DF,  p-value: < 2.2e-16
# Y superponemos la recta de ajuste
plot(x, logCCDF)
abline(fit$coefficients, lty=2, col="red")

Como siempre, la función summary da mucha información. Del ejemplo anterior, lo importante es saber identificar los parámetros ajustados (ordenada en el origen y pendiente). Localícelos sabiendo que en nuestro ejemplo la recta teórica, por simple inspección de \(F^C(x)\), debería tener ordenada en el origen nula y pendiente \(-\lambda\).

De toda esa información sobre el ajuste, también es importante la significancia de cada parámetro, es decir, cómo de buena es la regresión. La significancia se marca mediante el p-valor (a menor p-valor, mejor ajuste) y un código de asteriscos (a más asteriscos, mejor ajuste).

1.2 Propiedad sin memoria

Se dice que una variable aleatoria \(X\) no tiene memoria si cumple la siguiente propiedad:

\[\Pr(X > t+s \mid X>s) = \Pr(X > t)\]

Se puede confirmar que la variable aleatoria uniforme sí tiene memoria. Para ello, se va a comparar la función de supervivencia de la uniforme (parte derecha de la ecuación) con una función de supervivencia construida en dos pasos: (1) se seleccionan los valores mayores que un umbral dado \(s\), y (2) luego se les resta dicho valor.

Sea una uniforme entre 0 y 10, de la que obtenemos su función de distribución y, a partir de esta, la función de distribución complementaria (esto es, la función de supervivencia):

x <- runif(1000, 0, 10)
Fu <- ecdf(x)
plot(sort(x), 1-Fu(sort(x)))

Si se selecciona ahora el subconjunto de valores mayores de \(s=2\), se les resta dicha cantidad, y se repite el cálculo de la función complementaria, se puede apreciar la diferencia:

x2 <- x[x>2]-2
Fu2 <- ecdf(x2)
plot(sort(x), 1-Fu(sort(x)))
points(sort(x2), 1-Fu2(sort(x2)), col="red")

Si se repite el procedimiento pero para una variable aleatoria exponencial de media 5, se tiene que las funciones de supervivencia sí que se parecen mucho:

x <- rexp(1000, 1/5)
x2 <- x[x>5] - 5
Fu <- ecdf(x)
plot(sort(x), 1-Fu(sort(x)))
Fu2 <- ecdf(x2)
points(sort(x2), 1-Fu2(sort(x2)), col="red")

Lo que se puede confirmar con más rigor por medio del qqplot

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

1.3 Mínimo de variables exponenciales

Partimos de dos variables exponenciales \(X \sim \mathrm{Exp}(5)\) e \(Y \sim \mathrm{Exp}(10)\). En primer lugar, se crea la nueva variable \(Z \sim \min\{X, Y\}\). Para ello, se obtiene el mínimo de los dos vectores de muestras componente a componente con el comando pmin, como se hizo en la práctica anterior:

x <- rexp(1000, 1/5)
y <- rexp(1000, 1/10)
z <- pmin(x, y)

Se puede apreciar que las tres variables aleatorias ``parecen’’ seguir la función de densidad de una exponencial, donde la función es más plana cuanto mayor es la media (\(1/\lambda\))

plot(density(x), ylim=c(0,0.3))
lines(density(y), col='blue')
lines(density(z), col='red')

Se puede obtener la media de dicho mínimo (\(1/\lambda\)) con el comando mean

mean(z)
## [1] 3.252753

O bien se puede obtener la función de supervivencia en escala logarítmica, con lo que se logran dos objetivos: (1) confirmar que se trata de una línea (lo que supone un análisis más sólido que la estimación de la función de densidad) y (2) obtener el parámetro de la exponencial mediante una regresión lineal, ya que se trata de la pendiente de dicha recta (\(-\lambda\)):

Fx <- ecdf(z)
logCCDF <- log(1 - Fx(z))
# Se filtran valores "raros" (logaritmos de números negativos)
z <- z[is.finite(logCCDF)]
logCCDF <- logCCDF[is.finite(logCCDF)]

# Regresión lineal
fit <- lm(logCCDF ~ z)
plot(z, logCCDF)
abline(fit$coefficients, lty=2, col="red")

Y se comprueba que la pendiente coincide aproximadamente con el valor teórico: la nueva tasa es la suma de las tasas, \(\frac{1}{5}+\frac{1}{10}=0.3\).

-fit$coefficients["z"]
##         z 
## 0.3150273

2 Procesos de Poisson

2.1 Generación

Se va a comprobar que tanto la primera como la segunda definición de un proceso de Poisson sirven para generar dicho proceso. Para ello, en primer lugar se generará un proceso de llegadas siguiendo la segunda definición y se comprobará que cumple la primera definición, y en segundo lugar se realizará el camino inverso (se generará según la primera definición, y se comprobará que cumple la segunda definición).

2.1.1 2a definición -> 1a definición

Según la definición exponencial (segunda definición), un proceso de Poisson es aquel en el que el tiempo entre llegadas sigue una variable exponencial de media \(1/\lambda\). A continuación, partimos de una variable exponencial que representa el tiempo entre llegadas. La suma acumulativa produce un vector de tiempos absolutos de llegada:

inter <- rexp(1000, 2)
llegadas <- cumsum(inter)

Con este paso se termina la generación de un proceso de llegadas. Se puede comprobar que, dado que se han generado 1000 tiempos entre llegadas de media 1/2, la última llegada (la última posición del vector) se habrá producido cerca del instante 500

llegadas[length(llegadas)] # También serviría max(llegadas)
## [1] 524.2192

A continuación, se pasa a comprobar mediante la primera definición que dicho vector contiene instantes de llegada de un proceso de Poisson. Según esta definición, los eventos que se producen en un intervalo de longitud \(t\) siguen una variable discreta de Poisson de media \(\lambda t\).

Para obtener una secuencia con diferentes realizaciones del número de eventos de un intervalo de un tamaño dado, se realizan tres paso: 1) Se divide el tiempo total en intervalos de longitud 1 con el comando seq 2) Se asigna cada llegada a alguno de estos intervalos con el comando cut 3) Se cuenta el número de “asignaciones” por cada intervalo con tabulate

intervalos <- seq(0, max(llegadas), 1) # Intervalos de longitud 1
asignacion <- cut(llegadas, intervalos)
n_llegadas <- tabulate(asignacion)

Se puede comparar el histograma resultante con el valor teórico esperado mediante los siguientes comandos

hist(n_llegadas,probability=TRUE, breaks=seq(-0.5, 10.5, 1))
lines(0:10, dpois(0:10, 2), type = "h", lwd = 2, col="red")

También se podría comparar mediante un qqplot el vector resultante con uno generado mediante la función rpois (si bien el uso del qqplot con variables discretas no es habitual)

qqplot(n_llegadas, rpois(1000, 2 * 1))
abline(0, 1, lty=2, col="red")

Repita los comandos anteriores (generación de tiempos de llegadas y comparación con un proceso de Poisson) para diferentes valores de \(\lambda\) y \(t\).

2.1.2 1a definición -> 2a definición

Según la primera definición, dado un intervalo de tiempo de longitud \(t\) y una tasa de llegadas \(\lambda\), se puede emplear una variable de Poisson de media \(\lambda t\) para obtener un número de llegadas aleatorio en dicho intervalo.

N <- 1000 # Intervalos (no llegadas)
lambda <- 2
t <- 1
n <- rpois(N, lambda * t)

Para obtener en qué instante en concreto se ubican dichas llegadas, se sitúan al azar de forma uniforme (por el análisis de la distribución condicionada del instante de una llegada), que se ordenan de menor a mayor antes de incorporarlas al vector de llegadas general

llegadas <- c()
for (i in seq(1,N)) { 
  llegadas_intervalo <- runif(n[i], (i-1)*t, i*t)
  llegadas <- c(llegadas, sort(llegadas_intervalo))
}

A continuación, se pasa a comprobar mediante la segunda definición que dicho vector contiene instantes de llegada de un proceso de Poisson. Para ello, se calculan los tiempos entre llegadas y se comprueba que siguen una v.a. exponencial como la esperada

inter <- diff(llegadas)
mean(inter)
## [1] 0.486135
plot(density(inter), ylim=c(0,lambda))
curve(dexp(x, lambda), add=TRUE, lty=2, col="red")

Repita el proceso anterior para diferentes valores de \(\lambda\) y \(t\).

2.2 Agregación y descomposición

Por lo visto en el apartado anterior, a la hora de trabajar con procesos de Poisson suele ser más cómodo operar con los tiempos entre llegadas (esto es, las variables exponenciales). A continuación, se comprueban las propiedades de agregación y descomposición de procesos de Poisson, y el teorema de Palm-Khintchine.

2.2.1 Agregación

Se generan dos procesos de Poisson a diferentes tasas

N <- 2000
lambda1 <- 1
inter1 <- rexp(N, lambda1)
times1 <- cumsum(inter1)

lambda2 <- 3
inter2 <- rexp(N, lambda2)
times2 <- cumsum(inter2)

Dado que se ha generado el mismo número de muestras, un vector tiene una duración mayor que el otro en el tiempo, por lo que se limitan al mínimo:

max(times1); max(times2)
## [1] 2005.472
## [1] 668.2591
tmax = min(max(times1), max(times2))

times1 <- times1[times1 < tmax]
times2 <- times2[times2 < tmax]

Para agregar dos procesos, basta con unir los vectores y ordenarlos:

llegadas <- sort(c(times1, times2))

Se puede comprobar que los tiempos entre llegadas siguen una v.a. exponencial de tasa igual a la suma de las tasas

inter <- diff(llegadas)
mean(inter)
## [1] 0.2498687
plot(density(inter), ylim=c(0,lambda1+lambda2))
curve(dexp(x, lambda1+lambda2), add=TRUE, lty=2, col="red")

Y para más rigor, se puede comparar la distribución resultante con la esperada mediante un qqplot:

qqplot(inter, rexp(length(inter), lambda1 + lambda2))
abline(0, 1, lty=2, col="red")

2.2.2 Descomposición / Muestreo

Se genera un proceso de Poisson a tasa 4:

N <- 2000
lambda <- 4
inter <- rexp(N, lambda)
llegadas <- cumsum(inter)

Para seleccionar p.ej. 1/4 de las muestras al azar, primero se genera un vector de índices donde cada valor es cierto con probabilidad 1/4, y se seleccionan dichas llegadas:

prob <- 1/4
indices <- runif(length(llegadas)) < prob
llegadas_subset <- llegadas[indices]

Se puede comprobar que el tiempo entre llegadas sigue una v.a.exponencial de media la esperada, p.ej. con un qqplot:

interarrivals <- diff(llegadas_subset)
qqplot(interarrivals, rexp(length(interarrivals), lambda*prob))
abline(0, 1, lty=2, col="red")

2.2.3 Teorema de Palm-Khintchine

El tiempo entre llegadas en un proceso, en general, no tiene por qué ser exponencial. Sea por ejemplo un codificador de voz sobre IP (VoIP) que genera una trama cada 10~ms, pero con una variabilidad tal que el tiempo entre tramas sigue una variable aleatoria uniforme entre 8~ms y 12~ms.

interu <- runif(1000, 8, 12)
timesu <- cumsum(interu)

Obviamente la función de densidad de los tiempos entre llegada dista mucho de ser una variable aleatoria exponencial:

plot(density(interu))

E incluso si se cuentan el número de llegadas en ventanas de una longitud fija (por ejemplo, 15), la distribución del número de llegadas no sigue una variable de Poisson

intervalos <- seq(0, max(timesu), 15) # Intervalos de longitud 15
asignacion <- cut(timesu, intervalos)
n_llegadas <- tabulate(asignacion)
hist(n_llegadas,probability=TRUE, breaks=seq(-0.5, 4.5, 1))

Pruebe a cambiar la longitud del intervalo y compruebe el impacto en el histograma resultante

Sin embargo, el agregado de un número suficiente de procesos de llegada tiende a comportarse como un proceso de Poisson, por lo que el tiempo entre llegadas se parecerá a un tiempo entre llegadas exponencial. Sea el caso del agregado de cuatro procesos como el anterior

llegadasu1 <- cumsum(runif(1000, 8, 12))
llegadasu2 <- cumsum(runif(1000, 8, 12))
llegadasu3 <- cumsum(runif(1000, 8, 12))
llegadasu4 <- cumsum(runif(1000, 8, 12))

llegadas <- sort(c(llegadasu1, llegadasu2, llegadasu3, llegadasu4))

Se puede comprobar que el tiempo entre llegadas empieza a parecerse a una variable aleatoria exponencial, esto es, se “escora” hacia los valores más bajos.

interarrivals <- diff(llegadas)
plot(density(interarrivals), ylim=c(0,0.4))
curve(dexp(x, 0.4), add=TRUE, lty=2, col="red")

Si bien el qqplot muestra que para cuatro procesos el parecido con una variable aleatoria exponencial es relativamente bajo

qqplot(interarrivals, rexp(length(interarrivals), 0.4))
abline(0, 1, lty=2, col="red")

Se puede probar a agregar más procesos y comprobar el parecido con la correspondiente v.a. exponencial

Ntotal <- 10
RateUno <- 1/((8+12)/2) # Tasa de llegadas de un proceso
llegadas_total <- c()
for (i in seq(1,Ntotal)) { 
  llegadas <- cumsum(runif(1000, 8, 12))
  llegadas_total <- sort(c(llegadas_total, llegadas))
}


interarrivals <- diff(llegadas_total)
qqplot(interarrivals, rexp(length(interarrivals), Ntotal*RateUno))
abline(0, 1, lty=2, col="red")

Pruebe a cambiar el número de procesos agregados, así la variable aleatoria del tiempo entre llegadas (p.ej., los límites de la v.a. uniforme), para ver el impacto en la distribución de los tiempos entre llegadas

3 Resolución de problemas con R

3.1 Problema 2.7

El retardo en una red se distribuye según una variable aleatoria exponencial de media 10 ms. Calcule el retardo medio de los paquetes con un retardo mayor de 50 ms.

Basta con generar la variable aleatoria y calcular la media del subconjunto de valores mayores de 10

N <- 1000
lambda <- 1/50
x <- rexp(N, lambda)
x <- x[x>10]
mean(x)
## [1] 61.69479

3.2 Problema 2.9

El tiempo de vida de un disco duro (HD) se puede modelar con una variable aleatoria exponencial de media \(1/\lambda\). Se desea almacenar 2 TB de información y hay dos configuraciones para los discos duros: * En la configuración RAID~0 (Data Striping) no existe redundancia y la capacidad total es la suma de las capacidades: si cualquier disco se estropea, la información se pierde. * En la configuración RAID~1 (Mirroring) la misma información es copiada en todos los discos y la capacidad total coincide con la de un disco: sólo cuando ningún disco funciona, la información se pierde. Se precisa saber la probabilidad de que dicha información esté disponible tras un año de uso para los siguientes casos:

  1. 1 HD de 2TB, tiempo medio de vida \(1/\lambda =\) 2 años

Para este caso, basta con calcular el número de veces que una exponencial de media 2 años supera el año de vida

required_lifetime <- 1
lambda <- 1/2
N <- 10000
hd_1 <- rexp(N, lambda)
survived <- hd_1[hd_1 > required_lifetime]

length(survived)/N
## [1] 0.6081
  1. 2 HD de 1TB cada uno, RAID0, tiempo medio de vida \(1/\lambda =\) 1 año

En esta configuración se necesita que los dos discos duros hayan sobrevivido al menos un año

N <- 10000
required_lifetime <- 1
lambda <- 1
hd_1 <- rexp(N, lambda)
hd_2 <- rexp(N, lambda)

sobreviven = (hd_1 > required_lifetime) & (hd_2 > required_lifetime)
survived <- sum(sobreviven, na.rm = TRUE)
survived/N
## [1] 0.1312
  1. 2 HD de 2TB cada uno, RAID1, tiempo medio de vida \(1/\lambda = 1 year\)
N <- 10000
required_lifetime <- 1
lambda <- 1
hd_1 <- rexp(N, lambda)
hd_2 <- rexp(N, lambda)

sobreviven = (hd_1 > required_lifetime) | (hd_2 > required_lifetime)
survived <- sum(sobreviven, na.rm = TRUE)
survived/N
## [1] 0.6001

3.3 Problema 3.8

Una aplicación genera tramas según un proceso de Poisson a tasa \(\lambda=2\)~tramas/s. Para que el interfaz no entre en modo de ahorro de energía, un módulo de control genera una trama si no se detecta tráfico durante un tiempo \(Timeout=1\)~s. Calcule la tasa de transmisión de tramas, y discuta si se trata de un proceso de Poisson.

Una forma de simular el proceso de generación de tramas es partir de varias muestras de una variable aleatoria exponencial de media 1/2. Dichas muestras serán los tiempos entre llegadas si son menores de \(Timeout=1\), en caso contrario es preciso añadir una llegada. En primer lugar, se genera el proceso de llegadas de la aplicación

N <- 2000
lambda <- 2
inter <- rexp(N, lambda)
llegadas <- cumsum(inter)

A continuación, recorremos dicho proceso de llegadas, y en caso de que entre una llegada y la siguiente pase más de un tiempo \(Timeout\), se generan tantas tramas como instantes de tiempo quepan

Timeout <- 1
salidas_interfaz <- c()
t_now <- 0
for (t_llegada in llegadas) {
  if (t_llegada - t_now < Timeout) { # Llega antes del umbral
    salidas_interfaz <- c(salidas_interfaz, t_llegada)
  } else { # Llega tarde, hay que ver cuanto de tarde
    n_thresholds <- floor((t_llegada-t_now)/Timeout)
    t_generados <- t_now + seq(1,n_thresholds)*Timeout
    salidas_interfaz <- c(salidas_interfaz, t_generados, t_llegada)
  }
  t_now <- t_llegada
}

Se puede comprobar que el tiempo medio entre tramas es menor que en el caso anterior y la tasa de generación de tramas es mayor

iarrival_salida <- diff(salidas_interfaz)
mean(iarrival_salida)
## [1] 0.4235571
1/mean(iarrival_salida)
## [1] 2.360957

Se puede comprobar que el tiempo entre llegadas dista seguir una variable aleatoria exponencial, dado que hay una acumulación de valores en \(Timeout\).

plot(density(iarrival_salida), ylim=c(0, 1/mean(iarrival_salida)))
curve(dexp(x, 1/mean(iarrival_salida)), add=TRUE, lty=2, col="red")

Usando un qqplot la diferencia también es evidente:

qqplot(iarrival_salida, rexp(length(iarrival_salida), 1/mean(iarrival_salida)))
abline(0, 1, lty=2, col="red")