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”:

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

summary(x)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0014  1.6718  3.7766  5.3897  7.2926 42.3117

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

plot(density(x))

# f(x) teórica
curve(dexp(x, 1/5), 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/5), 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}\]

trabajamos con una exponencial, a la que podemos aplicar un logaritmo y obtener una recta cuya pendiente es la tasa de la exponencial:

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

La función de supervivencia o CCDF:

plot(x, 1-Fx(x))

Aplicando logaritmo, obtenemos una recta:

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

Y de aquí podemos obtener la tasa por regresión lineal. (¿Qué es una regresión lineal? Aquí un poco de teoría 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
# nos molesta al hacer la regresión y que debemos 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.1875699  0.0002165 -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

Partimos de una variable uniforme entre 0 y 10. Vamos a avanzar 2 unidades de tiempo, esto es, nos quedamos con las muestras mayores que 2 y restamos 2 a todas ellas:

x <- runif(1000, 0, 10)
x2 <- x[x>2] - 2

La nueva distribución x2 difiere de la original, luego la distribución uniforme tiene memoria:

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

Sin embargo, se puede demostrar que la distribución exponencial no tiene memoria:

x <- rexp(1000, 1/5)
x2 <- x[x>2] - 2
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), Y \sim \mathrm{Exp}(10)\). Buscamos la nueva variable \(Z \sim \min\{X, Y\}\). Para ello, necesitamos obtener el mínimo de los dos vectores de muestras componente a componente:

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

Si las variables anteriores conforman un componente único (por ejemplo un router), para estimar el parámetro \(\lambda\), ajustamos una regresión lineal sobre la función de supervivencia:

Fx <- ecdf(z)
logCCDF <- log(1 - Fx(z))
z <- z[is.finite(logCCDF)]
logCCDF <- logCCDF[is.finite(logCCDF)]

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

cuya 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

Poisson Process

Poisson Process

2.1 Generació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\). Es decir, dado un tiempo entre llegadas distribuido según una exponencial de tasa \(\lambda\), si contamos los eventos que suceden en cada intervalo \(t\), estos se distribuyen según una variable de Poisson de parámetro \(\lambda t\).

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)
times <- cumsum(inter)

El número de llegadas por unidad de tiempo estará distribuido según una variable de Poisson. Vamos a generar una función que permita contarlas para compararlas con una distribución \(\mathrm{Pois}(2)\):

# Divide un fragmento de tiempo en intervalos de tamaño "t" y cuenta el número de llegadas
count_arrivals <- function(times, t=1) {
  # número de intervalos de tamaño "t" (vector)
  breaks <- seq(0, max(times), t)

  # CUT divide el vector "times" en intervalos (especificados en "breaks")
  # TABULATE cuenta el número de veces que un intervalo ha tenido una llegada
  tabulate(cut(times, breaks))
}

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

Si duplicamos el intervalo de conteo, pasamos a una \(\mathrm{Pois}(4)\):

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

También podemos seguir el camino inverso. A partir de la primera definición, dado un intervalo de tiempo de longitud \(t\), podemos usar una variable de Poisson para determinar el número de llegadas correspondientes. Entonces, siguiendo la probabilidad condicionada de una llegada, simplemente podemos fijar aleatoriamente esas llegadas en dicho intervalo, y os tiempos entre llegadas resultantes seguirán una distribución exponencial.

lambda <- 2
t <- 1
x <- rpois(1000, lambda * t)

# Situamos de forma aleatoria llegadas x[i] en cada intervalo i
place <- function(i) {
  # generamos las llegadas
  arrivals <- runif(x[i], 0, t)

  # Almacenamos las llegadas en el intervalo correspondiente añadiendo i*tamaño del intervalos
  sort(arrivals) + (i-1)*t
}
# times format: 
#  - [[interval]] [1] arrival1 arrival2
times <- sapply(1:length(x), place)

# Extrae las llegadas de los intervalos en un vector, reemplazando el tiempo absoluto de la llegada por la diferencia entre llegadas.
inter <- diff(unlist(times))

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

2.2 Agregación y descomposición

En general, a la hora de trabajar con procesos de Poisson, es más cómodo partir de los tiempos entre llegadas (la variable exponencial subyacente). A continuación, partimos de dos procesos de Poisson. Comprobamos que ambos vectores están descompensados: el de menor tasa (mayor media) tiene una duración mayor en el tiempo, por lo que hay que recortar dicho vector. La agregación consiste simplemente en la concatenación y ordenación de los vectores de tiempo absolutos:

lambda1 <- 1
inter1 <- rexp(1000, lambda1)
times1 <- cumsum(inter1)

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

max(times1); max(times2)
## [1] 1002.422
## [1] 345.4732
times1 <- times1[times1 < max(times2)]
times2 <- times2[times2 < max(times1)]
times_agg <- sort(c(times1, times2))
inter_agg <- diff(times_agg)

Comprobamos que el proceso resultante sigue una distribución exponencial de tasa igual a la suma de tasas:

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

Ahora realizamos una descomposicón del proceso agregado de tasa 4. Para ello, cada llegada es descartada o no con probabilidad 0.5. El proceso resultante debería seguir una exponencial de tasa mitad de la original:

prob <- 1/2

# Descomposición usando un bucle for
# times_dec = c()
# for (i in 1:length(times_agg)){
#   random_val <- runif(1)
#   if (random_val < prob){
#      times_dec <- c(times_dec, times_agg[i])
#   }
# }

# Descomposición usando directamente vectores 
times_dec <- subset(times_agg, runif(length(times_agg)) < prob)

inter_dec <- diff(times_dec)
qqplot(inter_dec, rexp(length(inter_dec), (lambda1 + lambda2) * prob))
abline(0, 1, lty=2, col="red")

2.3 Teorema de Palm-Khintchine

Supongamos que un switch tiene conectados 64 teléfonos VoIP. Una llamada VoIP genera paquetes a tasa constante de forma que el tiempo entre paquetes sigue una distribución uniforme entre 8 y 12 ms. El proceso de llegadas de cada teléfono por separado no puede considerarse de Poisson, porque los tiempos entre llegadas no siguen una exponencial con la misma media:

inter <- runif(1000, 8, 12)
qqplot(inter, rexp(length(inter), 1 / ((8 + 12) / 2)))
abline(0, 1, lty=2, col="red")

Sin embargo, como asegura el teorema de Palm-Khintchine, la agregación de estas 64 fuentes de tráfico sí está bastante próxima a un proceso de Poisson:

f <- function(i){
  samples <- runif(1000, 8, 12)
  return(cumsum(samples))
}
# Matriz de 64 columnas y 1000 filas
traffic_data <- sapply(1:64, f)
vectorized_traffic_data <- unlist(traffic_data)

times_agg <- sort(vectorized_traffic_data)
inter_agg <- diff(times_agg)

qqplot(inter_agg, rexp(length(inter_agg), 64 / ((8 + 12) / 2)))
abline(0, 1, lty=2, col="red")

3 Resolución de problemas con R

3.1 Problema 1

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.

lambda <- 1/10
f <- function(x){
  x <- rexp(1000, lambda)
  return(mean(x[x>50]))
}
avg <- sapply(1:1000, f)
mean(avg)
## [1] 60.02601

3.2 Problema 2

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 se precisa saber la probabilidad de que dicha información esté disponible tras un año de uso, existiendo tres posibles configuraciones:

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

  1. 1 HD 500€, 2TB, \(1/\lambda = 2 years\)
expected_lifetime <- 1

f <- function(x){
  hd_1 <- rexp(10000, 1/2)
  greater <- hd_1[hd_1 > expected_lifetime]
  return(length(greater)/length(hd_1))
}

r <- sapply(1:1000, f)
print(mean(r))
## [1] 0.6065801
  1. 2 HD 100€, 1TB cada uno, \(1/\lambda = 1 year\), Raid 0
expected_lifetime <- 1

f <- function(x){
  hd_1 <- rexp(1000, 1)
  hd_2 <- rexp(1000, 1)
  greater <- 0
  for (i in 1:length(hd_1)){
    if((hd_1[i] > expected_lifetime) && (hd_2[i] > expected_lifetime)){
      greater <- greater + 1
    }
  }
  return(greater/length(hd_1))
}
r <- sapply(1:1000, f)
print(mean(r))
## [1] 0.135007
  1. 2 HD 200€, 2TB cada uno, \(1/\lambda = 1 year\), Raid 1
expected_lifetime <- 1

f <- function(x){
  hd_1 <- rexp(1000, 1)
  hd_2 <- rexp(1000, 1)
  greater <- 0
  for (i in 1:length(hd_1)){
    if((hd_1[i] > expected_lifetime) || (hd_2[i] > expected_lifetime)){
      greater <- greater + 1
    }
  }
  return(greater/length(hd_1))
}
r <- sapply(1:1000, f)
print(mean(r))
## [1] 0.600491

3.3 Problema 3

El tiempo entre dos peticiones de un usuario a YouTube se puede modelar con una variable aleatoria distribuida uniformemente entre 0 y 40 minutos. El 60% de las peticiones se realizan a servidores en EEUU. Sea un una celda de telefonía móvil sirviendo a 100 usuarios:

  • Caracterice el proceso de generación de peticiones de video a EEUU, indicando las suposiciones realizadas.

  • Calcule la probabilidad de que en un minuto se produzcan menos de dos peticiones a EEUU.

Este ejercicio está pensado para utilizar el teorema de Palm-Khintchine, para poder asumir que el proceso resultante es de Poisson, tiempo entre llegadas exponencial y tasa \(0.6\cdot 100 \cdot 2 / 40 = 3\) peticiones/minuto. Pero, ¿es realmente similar a un proceso de Poisson?, simularemos el proceso con la información que tenemos para intentar responder a la pregunta:

n <- 10000
users <- 100
t.min <- 0
t.max <- 40
prob <- 0.6
aggr_times <- function(x) {
    # Generamos un vector de muestras
    u <- cumsum(runif(n/users/prob, t.min, t.max))
    # seleccionamos aquellas muestras que van a servidores de EEUU
    return(u[runif(n/users/prob) < prob])
}
# Funcion para calcular tiempo entre llegadas
inter_times <- function(x) {
  aggr_times_data <- sapply(1:users, aggr_times)
  return(diff(sort(unlist(aggr_times_data))))
}

data <- inter_times()

Este proceso de llegadas parece ajustarse a una distribución exponencial:

# Comparamos con una exponendial
qqplot(data, rexp(length(data), 3))
abline(0, 1, lty=2, col="red")

Por tanto, podemos decir que las peticiones a EEUU son un muestreo del proceso, siendo entonces otro proceso de Poisson con tasa:

out <- sapply(1:1000, function(i) {
  # cell contiene los tiempos entre tramas
  cell <- inter_times()
  return(1/mean(cell))
})
t.test(out)
## 
##  One Sample t-test
## 
## data:  out
## t = 1712.2, df = 999, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  2.700414 2.706611
## sample estimates:
## mean of x 
##  2.703512

Por lo tanto, la probabilidad de tener menos de dos peticiones en un minuto:

out <- sapply(1:1000, function(i) {
  # count_arrivals was defined at Poisson processes section
  counts <- count_arrivals(cumsum(inter_times()), t=1)
  sum(counts < 2) / length(counts)
})
t.test(out)
## 
##  One Sample t-test
## 
## data:  out
## t = 629.01, df = 999, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.2682669 0.2699460
## sample estimates:
## mean of x 
## 0.2691065

3.4 Problema 4

Una aplicación genera tramas según un proceso de Poisson a tasa \(\lambda\)~tramas/s. Para que el interfaz inalámbrico 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 \(T=1/\lambda\)~s. Calcule la tasa de transmisión de tramas, y discuta si se trata de un proceso de Poisson.

En primer lugar, implementaremos una función que nos permita representar ese proceso. Para ello, comenzamos generando muestras distribuidas exponendialmente que representarán el tráfico que proviene de la aplicación. Tras esto, pondremos tramas de control en aquellos huecos de tamaño \(T\ge1/\lambda\):

rgen <- function(n, lambda) {
  # samples generation
  # generación de muestras
  e1 <- rexp(n, lambda)
  # Buscamos los huecos con tamaño T>=1/lambda
  gap <- which(e1 >= 1/lambda)
  e1 <- as.list(e1)
  # Obtenemos el número de tramas de control que deberíamos generar
  for (j in gap) {
    N <- as.integer(e1[[j]] * lambda)
    reps <- rep(1/lambda, N)
    e1[[j]] <- c(reps, e1[[j]] - sum(reps))
  }
  unlist(e1)[1:n]
}

Finalmente, obtenemos la tasa para varias simulaciones. Asumimos, por ejemplo, \(\lambda=2\):

out <- sapply(1:1000, function(i) {
  1/mean(rgen(1000, 2))
})
t.test(out)
## 
##  One Sample t-test
## 
## data:  out
## t = 1705.2, df = 999, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  3.163797 3.171087
## sample estimates:
## mean of x 
##  3.167442

Es evidente que el proceso resultante no es de Poisson. En la función de densidad podemos ver un pico en la parte derecha que explica el tráfico que proviene del módulo de control.

e1 <- rgen(1000, 2)
plot(density(e1))
curve(dexp(x, 2), add=TRUE, lty=2, col="red")

Usando qqplot para compararlo con una distribución exponencial, vemos lo mismo:

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