Introducción

En este LAB simularemos y veremos algunas propiedades de los procesos de Poisson vistos en el teórico.

Con el abordaje de los Procesos de Poisson estudiamos la ocurrencia de eventos en el tiempo o en el espacio. En el caso de variables continuas, fijamos un valor umbral y modelamos el número de veces que se supera dicho valor en un período de tiempo dado (High level excedance).

Recordemos que en la primera parte del curso, mediante el enfoque de Fisher-Tippett-Gnedenko basado en las distribuciones extremales, nos ocupamos de modelar los valores máximos de cierto fenómeno, es decir, para un período de tiempo dado, considerábamos el valor máximo de variables aleatorias \(iid\). Con los Procesos de Poisson (PP) no nos enfocamos en “lo peor” o lo más extremo, sino, en eventos extremos o inusuales que superan cierto umbral definio por el/la investigador/a. Además nos permite trabajar con la ocurrencia de eventos (discretos).

Ejemplos :

¿Otros ejemplos?

Comenzaremos realizando un repaso de las distribuciones de Poisson y Exponencial ya que son fundamentales para entender los Procesos de Poisson. Luego, continuaremos trabajando con los procesos de Poisson en este orden:

1- Procesos de Poisson simples (PP),

2- Procesos de Poisson Compuestos (PPC),

3- Procesos de Poisson No Homogéneos (PPNH).

Distribucion de Poisson

La distribución de Poisson es una distribución de probabilidad discreta que expresa la probabilidad de que ocurra un determinado número de eventos durante cierto período de tiempo.

Sea \(\lambda > 0\) y \(X\) una variable aleatoria discreta. Si \(X\) tiene una distribución de Poisson con parámetro \(\lambda\), entonces notaremos \(X \sim {\sf Poisson}(\lambda)\).

Entonces, si \(X \sim {\sf Poisson}(\lambda)\) su función de probabilidad es

\[\begin{equation} P[X=k] = \frac{\ e^{-\lambda}\lambda^{k}}{\ k!} \end{equation}\]

donde \(k=0,1,2,...\) es el número de ocurrencias del evento o fenómeno.

Nota: Recordemos que la en la distribución de Poisson, \(X\) es una variable aleatoria discreta, es decir, sólo puede tomar valores enteros \(k=0,1,2,...\), por lo que tiene sentido hablar de la probabilidad de que \(X\) sea exactamente igual a \(k\) y en vez de expresar \(P[X\leq k]\) que corresponde a la función de distribución acumulada (FDA) para el caso de que \(X\) sea una v.a. continua.

Por su parte, el parámetro \(\lambda\) representa el número de veces que se espera que ocurra el fenómeno durante un intervalo dado. Tanto el Valor Esperado como la Varianza de una variable aleatoria \(X\) con distribución de Poisson son

Entonces \(\lambda\) representa la tasa promedio de ocurrencias del fenómeno en el intervalo considerado, con ua desviación estándar de \(\sqrt{\lambda}\).

Simulación

Simulemos en R la función de probabilidad de una variable aleatoria tal que \(\;X \sim {\sf Poisson}(5)\). La podemos expresar como

\[ P[X = k]=e^{-5}\frac{5^k}{k!},\quad k=0,1,2,... \] para cada valor de \(k\) se calcula la probabilidad.

lmbd5=5 #parámetro lambda
x=0:50 # son los posibles valores de la variable X (sus cuantiles o puntos de soporte)

En esta simulación \(X=1,\dots, 50\). Elegimos cortar el código en 50 por que no podemos trabajar con valores infinitos, entonces elegimos un rango de valores lo suficientemente grande. Es decir, cortar en 50 nos basta para graficar casi toda la masa de probabilidad.

p5=dpois(x,lmbd5) # probabilidad
p5
##  [1] 6.737947e-03 3.368973e-02 8.422434e-02 1.403739e-01 1.754674e-01
##  [6] 1.754674e-01 1.462228e-01 1.044449e-01 6.527804e-02 3.626558e-02
## [11] 1.813279e-02 8.242177e-03 3.434240e-03 1.320862e-03 4.717363e-04
## [16] 1.572454e-04 4.913920e-05 1.445271e-05 4.014640e-06 1.056484e-06
## [21] 2.641211e-07 6.288597e-08 1.429227e-08 3.107014e-09 6.472947e-10
## [26] 1.294589e-10 2.489595e-11 4.610361e-12 8.232787e-13 1.419446e-13
## [31] 2.365743e-14 3.815715e-15 5.962055e-16 9.033417e-17 1.328444e-17
## [36] 1.897777e-18 2.635801e-19 3.561893e-20 4.686701e-21 6.008592e-22
## [41] 7.510739e-23 9.159438e-24 1.090409e-24 1.267918e-25 1.440816e-26
## [46] 1.600906e-27 1.740116e-28 1.851187e-29 1.928320e-30 1.967673e-31
## [51] 1.967673e-32

Nos devuelve los valores de \(P[X = 0],\dots,P[X = 50]\).

plot(p5, type = "h", lwd = 2,
     main = "Función de probabilidad", ylab="P(X=x)",xlab="Número de eventos")

sum(p5)
## [1] 1

Eso confirma que cortar en 50 fue una excelente elección para representar toda la masa de probabilidad.

Veamos ahora al comportamiento de la función de probabilidad considerando diferentes lambdas:

lmbd10 = 10
lmbd20 = 20

p10 = dpois(x, lmbd10)
p20 = dpois(x, lmbd20)
plot(x, p5,  type="h", lwd=2, col="black",
     main="Función de probabilidad Poisson(λ)",
     ylab="P(X = x)", xlab="Número de eventos")
lines(x, p10, type="h", lwd=2, col=rgb(1,0,0,0.7))
lines(x, p20, type="h", lwd=2, col=rgb(0,1,0,0.7))
legend("topright", legend=c("λ=5", "λ=10", "λ=20"),
       lty=1, col=c("black", rgb(1,0,0,0.7), rgb(0,1,0,0.7)),
       lwd=2, box.lty=0)

Por tanto, al aumentar \(\lambda\):

  • La distribución se desplaza hacia la derecha (la media aumenta).
  • La dispersión también aumenta (más valores posibles con probabilidad no nula).
  • La forma se va volviendo más simétrica (por el Teorema Central del Límite, Poisson \(\approx\) Normal cuando \(\lambda\) grande).

Considerando una \(X \sim {\sf Poisson}(5)\), ¿cuál es la probabilidad de observar más de cinco eventos en un intervalo de tiempo? Es decir, ahora calculamos probabilidades del tipo \(P[X>k]\), con \(k=5\), es decir, la probabilidad de que ocurran más de 5 eventos en ese intervalo.

En este esquema el cálculo teórico de la probabilidad \(P[X > 5]\) es

\[ P[X > 5] = 1 - P[X \le 5]. \]

La cantidad \(P[X \le 5]\) corresponde a la función de distribución acumulada (CDF) evaluada en 5, es decir

\[ P[X \le 5] = \sum_{k = 0}^{5} e^{-\lambda} \frac{\lambda^k}{k!}. \]

Por lo tanto, si \(X \sim \text{Poisson}(\lambda = 5)\), la probabilidad de que ocurran más de cinco eventos en un intervalo de tiempo es

\[ P[X > 5] = 1 - \sum_{k = 0}^{5} e^{-5} \frac{5^k}{k!}. \]

# distribution function
ppois(5,lambda=5,lower.tail=FALSE) # lower.tail = FALSE porque P(X>q)
## [1] 0.3840393

¿Qué valor toma la probabilidad de observar más de cinco eventos en un intervalo de tiempo si consideramos que la variable \(X \sim {\sf Poisson}(10)\)?

ppois(5,lambda=10,lower.tail=FALSE)
## [1] 0.932914

Es decir, hay aproximadamente un 38.4% de probabilidad de observar más de cinco eventos en un intervalo, cuando en promedio esperamos 5.

Ahora, ¿cuál la probabilidad de obtener exactamente 5 eventos (\(P[X=5]\))? Consideramos una \(X \sim {\sf Poisson}(5)\)

p_menori4 = ppois(4, lambda = 5)                   # P(X ≤ 4)
p_mayor5  = ppois(5, lambda = 5, lower.tail = FALSE) # P(X > 5)
1 - p_menori4 - p_mayor5                           # 1 - [P(X ≤ 4) + P(X > 5)]
## [1] 0.1754674

La función de distribución acumulada (CDF) se define como:

\[ P[X \le k] = \sum_{i = 0}^{k} P[X = i]. \]

Por lo tanto, la probabilidad puntual de obtener exactamente cinco eventos se puede expresar como la diferencia entre dos valores acumulados consecutivos como

\[ P[X = 5] = P[X \le 5] - P[X \le 4]. \]

Además, como se cumple que

\[ P[X \le 4] + P[X > 5] + P[X = 5] = 1, \]

podemos despejar \(P[X = 5]\) y obtener una expresión equivalente:

\[ P[X = 5] = 1 - \big[ P[X \le 4] + P[X > 5]\big]. \] Esta es una forma indirecta pero correcta de aislar la probabilidad de \(X = 5\), utilizando únicamente la función de distribución acumulada \(P[X \le k]\).

Alternativamente, obtenemos el mismo resultado con

dpois(5, lambda = 5)
## [1] 0.1754674

Nos devuelve la función de probabilidad puntual \(P[X=5]\).

La siguiente función es útil para simular variables aleatórias con distribución de Poisson.

set.seed(71)
rpois(70,lambda=3) #genero 70 números aleatorios de una variable X con distribución de Poisson de parámetro lambda=3
##  [1] 2 3 2 2 2 6 4 5 2 3 1 4 1 4 5 4 1 4 5 6 3 2 4 2 1 3 2 2 2 5 1 3 4 1 4 5 3 2
## [39] 3 0 2 3 0 4 2 3 3 1 0 5 0 1 3 3 6 8 3 3 4 1 4 2 1 1 2 1 3 5 3 0

Dicha variable podría representar por ejemplo:

  • Número de accidentes de tránsito en períodos de una semana

  • Número de pasajes de autos en un cruce de calles determinado por minuto

  • Número de cirugías de urgnecia en cierto centro de salud por semana

  • Número de rachas de viento de más de 10 s de duración que superen 80 km/h en período de un mes, etc.

Distribucion Exponencial

La distribución exponencial es una distribución continua que se utiliza para modelar tiempos de espera para la ocurrencia de un cierto evento. Dicho de otro modo, es una distribución que modela el tiempo (o espacio) entre dos eventos de un proceso de Poisson cuando los eventos ocurren de manera independiente y a una tasa constante.

El parámetro \(\lambda\) representa la tasa promedio de ocurrencia de los eventos (número esperado de eventos por unidad de tiempo).

Una variable aleatoria continua tiene una distribución exponencial con parámetro \(\lambda > 0\) (\(X \sim {\sf Exp}(\lambda)\)) si su función de densidad es

\[\begin{equation} f_x (x) = \lambda e^{-\lambda x} \text{para\;} x \ge 0. \end{equation}\]

El Valor Esperado de una variable aleatoria continua \(X\) con distribución de Exponencial es \(E(X)=1/\lambda\) y la Varianza \(Var(X)=1/\lambda^2\).

Por lo tanto, el parámetro \(\lambda\) controla tanto la frecuencia de los eventos como la dispersión de los tiempos de espera:

Simulación

Simulemos y grafiquemos en R la función de densidad de una variable aleatoria \(X\) con distribución Exponencial considerando distintos \(\lambda=\{ 0.5, 1,2,10\}\).

x <- seq(0, 8, 0.1)
x
##  [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8
## [20] 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7
## [39] 3.8 3.9 4.0 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5.0 5.1 5.2 5.3 5.4 5.5 5.6
## [58] 5.7 5.8 5.9 6.0 6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 7.0 7.1 7.2 7.3 7.4 7.5
## [77] 7.6 7.7 7.8 7.9 8.0
exp0_5=dexp(x, rate =0.5)
exp1=dexp(x, rate = 1)
exp2=dexp(x, rate = 2)
exp10=dexp(x, rate = 10)
# lambda = 0.5
plot(x, exp0_5, type = "l",
     ylab = "", lwd = 2, col = "red")
# lambda = 1
lines(x, exp1, 
     ylab = "", lwd = 2, col = "orange")
# lambda = 2
lines(x, exp2, col = "blue", lty = 1, lwd = 2)
# lambda = 10
lines(x, exp10, col = "green", lty = 1, lwd = 2)
# Leyenda
legend("topright", c(expression(paste(lambda)),"0.5","1","2","10"),
       lty = c(0, 1, 1,1,1), col = c(NULL, "green","red","orange","blue"), box.lty = 0, lwd = 2)

  • Cuando \(\lambda\) = 0.5 (roja): La curva disminuye lentamente y los tiempos entre eventos son largos. Hay alta probabilidad de que el evento ocurra mucho después.

  • Cuando \(\lambda\) = 1 (naranja): La curva cae más rápido y los eventos son más frecuentes, el tiempo medio entre ellos es 1 unidad.

  • Cuando \(\lambda\) = 2 (azul): Cae aún más rápido y los eventos se concentran en tiempos pequeños.

  • Cuando \(\lambda\) = 10 (verde): La densidad está muy pegada al eje vertical y casi todos los eventos ocurren muy rápidamente después del anterior.

set.seed(71)
rexp(70,rate=3) #genero 70 números aleatorios de una variable X con distribución Exponencial de parámetro lambda=3
##  [1] 0.34161981 0.03673591 0.33420901 0.53513766 0.10780931 0.07809687
##  [7] 0.72590899 0.17793893 0.22809202 0.18449957 0.48563854 0.20395740
## [13] 0.14214603 0.28103122 0.13883231 0.50670574 0.45335732 0.33881018
## [19] 0.29643503 0.01212275 0.19761962 1.02970164 0.47119800 0.39963671
## [25] 0.07961654 0.98596914 0.10640556 0.43789230 0.09225297 0.03134409
## [31] 0.69863177 0.46297981 0.06909505 0.00266523 0.10998945 0.08472040
## [37] 0.10514211 0.62857246 0.10375860 0.47428787 0.33857677 0.61654458
## [43] 0.01855822 0.21054388 0.23962139 0.13853652 0.14485412 0.05467695
## [49] 0.06864550 0.20505334 0.46528114 0.55786762 0.32665841 0.82013817
## [55] 0.09163054 0.54620434 0.76922357 0.25175766 0.08412288 0.44317290
## [61] 0.41080195 0.33185927 0.80157849 0.13601630 0.25729766 0.13092189
## [67] 0.22726978 0.24763711 0.24002005 0.25742620

Relación entre los Procesos de Poisson y la Exponencial

Recordemos qué es un Proceso de Poisson

Un proceso de Poisson con tasa \(\lambda >0\) (denotado \({\sf PP}(\lambda)\)) es un proceso de conteo \(\{N(t), t \ge 0\}\) que cumple:

  1. \(N(0) = 0\)
  2. Incrementos independientes: los números de eventos en intervalos disjuntos son independientes.
  3. Incrementos estacionarios: el número de eventos en cualquier intervalo de longitud \(t\) tiene distribución Poisson:

\[ P[N(t+s) - N(s) = k] = e^{-\lambda t} \frac{(\lambda t)^k}{k!}. \]

De esto se deduce que el número esperado de eventos en un intervalo de longitud \(t\) es

\[ \mathbb{E}[N(t)] = \lambda t. \]

Definimos los tiempos entre eventos

Sean los tiempos de ocurrencia de los eventos

\[ 0 < T_1 < T_2 < T_3 < \dots \]

y los tiempos intereventos

\[ X_i = T_i - T_{i-1}, \quad i = 1, 2, 3, \dots \] con \(T_0 = 0\).

Estos \(X_i\) representan los tiempos de espera entre eventos sucesivos.

Distribución de los tiempos intereventos

Una propiedad fundamental del proceso de Poisson es que

\[ X_1, X_2, X_3, \dots \stackrel{\text{iid}}{\sim} {\sf Exponencial}(\lambda) \]

Es decir

  • Los tiempos entre eventos son iid.
  • Cada uno tiene distribución exponencial con el mismo parámetro \(\lambda\).

Por tanto, su función de densidad es \(f_X(x) = \lambda e^{-\lambda x}, \quad x \ge 0\).

Entonces

  • Si el número de eventos en un intervalo de longitud \(t\) sigue una distribución Poisson\((\lambda t)\),
    entonces el tiempo entre eventos sigue una distribución Exponencial\((\lambda)\).

  • Cuanto mayor es \(\lambda\), más frecuentes son los eventos, es decir, los tiempos intereventos más cortos.

  • Cuanto menor es \(\lambda\), los eventos son más raros y tiempos de espera más largos.

Procesos de Poisson

De modo general, podemos decir que un Proceso de Poisson es un proceso estocástico que registra la ocurrencia de sucesos en el tiempo.

Proceso de Poisson simple (PP(\(\lambda\)))

Si \(N\) es un proceso de conteo y \(\lambda\) \(>\) \(0\) diremos que \(N\) es un Proceso de Poisson de parámetro \(\lambda\) (\(N\) es PP(\(\lambda\))) si se cumple:

  1. Para todo intervalo \(J\) de los reales positivos, \(N(J)\) es una variable aleatoria que tiene distribución de Poisson de parámetro \(\lambda *longitud(J)\).

Esto significa que si tomamos cualquier intervalo de tiempo, por ejemplo \(J = [a,b]\),
el número de eventos que ocurren en ese intervalo es

\[ N(J) = N(b) - N(a), \]

y sigue una distribución de Poisson tal que

\[ N(J) \sim {\sf Poisson}(\lambda \cdot (b - a)). \]

El parámetro \(\lambda > 0\) es la tasa promedio de ocurrencia,
es decir, el número esperado de eventos por unidad de tiempo.

En consecuencia, el número esperado de eventos en el intervalo \(J\) es

\[ \mathbb{E}[N(J)] = \lambda \, \text{longitud}(J). \]

Entonces, el proceso genera eventos de manera constante en promedio. Por ejemplo, si se duplica la longitud del intervalo, se espera aproximadamente el doble de eventos.

  1. Si \(J\), \(L\), \(M\)….es una cantidad arbitraria de intervalos de reales positivos DISJUNTOS entonces \(N(J)\), \(N(L)\), \(N(M)\) son variables aleatorias independientes.

Esto quiere decir que si consideramos distintos períodos de tiempo sin superposición (por ejemplo, la mañana y la tarde), el número de eventos ocurridos en esos períodos no se afectan entre sí. Es decir, el número de eventos entre 0 y 1 hora no da información sobre cuántos habrá entre 1 y 2 horas. El proceso no tiene memoria: los conteos en distintos intervalos son independientes.

Simulación

Simulemos un Proceso de Poisson simple. Vamos a generar una trayectoria simulada del proceso de Poisson \(N(t)\), es decir, una función escalonada que muestra cómo aumenta el número de eventos en el tiempo.

k = 50        # número total de eventos que vamos a simular
lambda = 2    # tasa del proceso de Poisson
set.seed(15)  # para reproducibilidad

Ahora, generamos los de tiempos entre eventos: \(X_i\sim \text{Exp}(\lambda)\).

expon=rexp(k,rate=lambda) #genero 50 números aleatorios de una variable X con distribución Exponencial de parámetro lambda=2

Ahora, obtenemos los tiempos de ocurrencia. Se calculan las sumas acumuladas \(T_i=X_1+\dots+X_i\). Cada \(T_i\) representa el tiempo exacto en que ocurrió el evento número \(i\).

Tiempo = cumsum(expon) #calculo las sumas acumuladas de dichas observaciones lo que se va a corresponder con los tiempos de ocurrencia de los eventos

Definimos \(N(t)=0,1,...,k\) y graficamos

N = 0:k #proceso de conteo
#grafico el tiempo en función del número de eventos 
plot(stepfun(Tiempo,N),do.points=F,xlab="t",ylab="N(t)")

  • Cada salto vertical representa un evento nuevo.

  • Cada tramo horizontal representa un período sin eventos.

  • El resultado es una típica trayectoria del proceso de Poisson.

Si bien aquí ya sabemos cual es nuestro \(\lambda\) pues nuestros datos son generados a partir de una simulación en la cual definimos \(\lambda=2\), hagamos el ejercicio de estimar lambda (\(\hat \lambda\)) manualmente a partir de los datos simulados.

Recordemos que los tiempos interevento tienen media \(1/\lambda\) ya que esta variable se distribuye según una exponencial de parámetro \(\lambda\).

Debemos en primer lugar, a partir de la variable Tiempo, calcular los tiempos intereventos, esto lo hacemos restando a cada tiempo el tiempo previo. En R esto lo hacemos fácilmente con la función diff().

Dif_tiempo=diff(Tiempo)
Dif_tiempo
##  [1] 0.973322890 0.127217424 0.088018729 0.331429021 1.377075686 0.146150910
##  [7] 0.009090392 0.206628568 0.155072176 0.048931711 0.303727399 0.293345952
## [13] 0.561699827 1.582367033 1.091146698 0.435553825 0.387287635 1.144275468
## [19] 0.158533380 0.364811196 0.534129738 0.461039446 0.343898076 2.275275105
## [25] 0.166700071 0.214119947 0.268503796 0.085963634 0.127862090 0.276130319
## [31] 0.317804436 0.393624199 0.490166232 0.173936632 0.274178611 0.508408451
## [37] 0.333017880 0.235766450 0.656200084 0.180911102 0.090463172 0.725742209
## [43] 0.141289949 0.264990433 0.056251295 0.581560194 0.298936418 0.004207106
## [49] 1.271231635

A partir de las diferencias, calculamos la media y computamos su inverso para obtener \(\lambda\).

mean(Dif_tiempo)
## [1] 0.4395509
lambda1=1/mean(Dif_tiempo)
lambda1
## [1] 2.275049

Observamos que el valor calculado de \(\lambda\) se encuentra bastante próximo a dos (valor teórico), a medida que aumentamos el número de datos (eventos) se aproximan cada vez más a este valor.

Probémoslo aumentando el número de eventos que registro. Tengamos en cuenta que en la vida real, para aumentar el número de registros deberíamos incrementar el esfuerzo de muestreo (mayor tiempo registrando el proceso).

Vamos a generar \(k_2=100000\) tiempos entre eventos con distribución Exponencia con distribución exponencial y \(\lambda=2\). Cada observación representa el tiempo que transcurre entre dos eventos sucesivos.

Si \(\lamda=2\) entonces \(E[X]=1/\lambda=1/2=0.5\). Entonces, en promedio, cada evento ocurre 0.5 unidades de tiempo después del anterior.

k2 = 100000  # defino k, número de eventos 
set.seed(122) 
expon2=rexp(k2,rate=lambda) 
expon2[0:5]
## [1] 0.06604132 1.54089157 0.07907347 0.12571680 0.17937594

Calculamos de los tiempos de ocurrencia acumulados

Tiempo2 = cumsum(expon2) #calculo las sumas acumuladas de dichas observaciones lo que se va a corresponder con los tiempos de ocurrencia de los eventos
Tiempo2[0:5]
## [1] 0.06604132 1.60693289 1.68600635 1.81172315 1.99109909
Dif_tiempo2=diff(Tiempo2)

Vuelve a darte los tiempos entre eventos, solo que pierde el primer valor (porque no hay \(T_0\)).

Ahora, calculamos \(\hat{\lambda}=1/\bar{X}\)

media=mean(Dif_tiempo2)
lambda2=1/media
lambda2
## [1] 1.998087

Son dos formas equivalentes de estimar la tasa \(\lambda\) de un proceso de Poisson homogéneo:

  1. A partir de los tiempos entre eventos, donde \(X_i \sim {\sf Exp}(\lambda)\):

\[ \hat{\lambda}_1 = \frac{1}{\bar{X}} \]

  1. A partir del número total de eventos observados \(N(0,t)\) y el tiempo total de observación \(t\):

\[ \hat{\lambda}_2 = \frac{N(0,t)}{t} \]

Cuando el número de eventos simulados es grande (\(k_2 = 100{,}000\)),
ambos métodos producen estimaciones muy cercanas al verdadero valor de \(\lambda\),
en concordancia con la Ley de los Grandes Números.

Gracias a la ley de los grandes números tenemos otra forma de estimar \(\lambda\) cuando \(t\) es grande:

\[N(0,t)/t\]

Es decir, el número total de eventos observados en el período \((0,t)\) sobre el tiempo total de observación.

k2/max(Tiempo2)
## [1] 1.998104

Nuevamente, asumiendo que desconocemos la estructura de los datos simulados, apliquemos el test de Lilliefors para determinar si los tiempos interevento siguen una distribución Exponencial. Recordemos que esta es una condición necesaria para asumir que estamos frente a un Proceso de Poisson simple (PP).

Usaremos la función LcKS del paquete {KScorrect}:

##install.packages("KScorrect")
require(KScorrect)
TestLillie=LcKS(Dif_tiempo, "pexp", nreps = 4999)
TestLillie$p.value
## [1] 0.3584

No rechazo \(H_0\) por lo que puedo asumir razonablemente que mis tiempos intereventos siguen una distribución Exponencial (cosa que en realidad ya sabíamos).

Proceso de Poisson Compuesto (PPC (\(\lambda,~G\)))

El Proceso de Poisson Compuesto (PPC) surge a partir de la necesidad de modelar fenómenos en donde los eventos no ocurren necesariamente uno a la vez, si no que pueden ocurrir más de un evento simultáneamente.

Repasemos la definición de este tipo de procesoÑ

Un Proceso de Poisson Compuesto (PPC) combina dos ideas fundamentales:

  1. Los tiempos de ocurrencia de los eventos siguen un proceso de Poisson de tasa \(\lambda\).
  2. La magnitud de cada evento (es decir, cuántas “ocurrencias” se producen en ese momento) es una variable aleatoria discreta \(S_i\) con distribución \(G\).

Entonces, definimos el proceso compuesto como

\[ M(t) = \sum_{i = 1}^{N(t)} S_i, \]

donde - \(N(t)\) es el número de eventos del proceso base hasta el tiempo \(t\), - y \(S_i\) representa cuántas ocurrencias simultáneas tiene el evento \(i\).

De este modo, mientras que el proceso de Poisson simple aumenta de uno en uno,
el proceso compuesto puede incrementarse en saltos de tamaño aleatorio.

Si \(N\) es un Proceso de Poisson de parámetro \(\lambda\) \(>\) \(0\), \(G\) es una distribución de probabilidad en los naturales (1,2,3,…), consideramos \(S_1,S_2,S_n,...iid\) con distribución \(G\) y construímos un nuevo proceso de conteo \(M\) de la forma siguiente:

Cuando \(N\) tiene su primer evento, \(M\) tiene \(S_1\) eventos simultáneos; Cuando \(N\) tiene su segundo evento, \(M\) tiene \(S_2\) eventos simultáneos…..(y así sucesivamente)

Decimos que M es un Proceso de Poisson Compuesto de parámetro \(\lambda\) \(>\) \(0\) y distribución de eventos \(G\) (\(M\) es \(PPC( \lambda,~G)\)).

En este tipo de proceso el tiempo medio interevento sigue siendo \(1/\lambda\), pero la tasa de incidencia media de eventos ahora es \(E(X)=\lambda ~E(G)\)

Simulación

Consideremos ahora los tiempos del primer ejemplo de PP y adjudiquemos a cada tiempo un número de eventos. Para esto último, a fin de simplificar, simularemos una distribución de Poisson y le sumaremos uno para evitar los ceros (es decir, garantizar al menos 1 evento por llegada).

# genera 50 números aleatorios que siguen una distribución de Poisson con media 1
nro_eventos=rpois(50,1)+1
nro_eventos #contiene los tamaños de los grupos o magnitudes de los 50 eventos del proceso base
##  [1] 2 1 3 2 1 1 1 1 1 1 3 2 2 2 2 2 2 2 4 3 1 1 3 3 2 2 1 1 1 1 1 2 2 2 2 4 2 1
## [39] 4 3 2 1 3 1 2 1 1 2 1 2
PPC=cbind(Tiempo,nro_eventos)
colnames(PPC)=c("Tiempo","Orden")
head(PPC,10)
##         Tiempo Orden
##  [1,] 0.102114     2
##  [2,] 1.075437     1
##  [3,] 1.202654     3
##  [4,] 1.290673     2
##  [5,] 1.622102     1
##  [6,] 2.999178     1
##  [7,] 3.145329     1
##  [8,] 3.154419     1
##  [9,] 3.361048     1
## [10,] 3.516120     1
plot(PPC[,"Tiempo"],cumsum(PPC[,"Orden"]), type="s")

Observamos que la gráfica muestras escalones o “saltos” que no son siempre unitarios, lo que nos debería inducir a pensar que estamos frente a un PPC.

Luego estimamos de forma empírica la función \(G\) como las frecuencias relativas de los Ordenes.

# Estimo la función G a partir de las frecuencias relativas.
G<- table(PPC[,"Orden"])/sum(table(PPC[,"Orden"]))
G
## 
##    1    2    3    4 
## 0.40 0.40 0.14 0.06

Podemos decir que estamos ante un \(PPC(\lambda=2, ~G)\).

Podríamos suponer que estos datos corresponden a los tiempos en que cierta variable supera determinado umbral, y en los casos en que este umbral es superado en “cluster” por más de un evento, se registra el número de eventos que ocurren juntos (Orden).

En este ejemplo, si nos interesa saber cual es el número de eventos promedio de Orden 4 por unidad de tiempo, debemos utilizar la estimación de \(\lambda= 2\) multiplicado por la probabilidad de de tener un evento de Orden 4 (\(G(4)=0.06\)) que nos queda \(0.06*2= 0.12\).

Los tiempos intereventos los tenemos calculados de la sección de PP (Dif_Tiempo) así como su promedio y a partir de allí el valor de lambda (que en realidad sabíamos previamente que tenía un valor porque lo simulamos los datos). Sabemos a su vez que el test de Lilliefors no rechaza la distribución Exponencial.

Algo que no hicimos previamente y que es importante, es evaluarlos tiempos intereventos mediante test de aleatoriedad para determinar que los datos sean iid. Para realizar los test de aleatoriedad utilizaremos el paquete {randtest} que tiene una batería de pruebas de aleatoriedad. Porbaremos varias. Si todas retienen H0 nos quedamos tranquilos de que los datos son iid.

##install.packages("randtests")
require(randtests)
runs.test(Dif_tiempo)
## 
##  Runs Test
## 
## data:  Dif_tiempo
## statistic = -1.1672, runs = 21, n1 = 24, n2 = 24, n = 48, p-value =
## 0.2431
## alternative hypothesis: nonrandomness
cox.stuart.test(Dif_tiempo)
## 
##  Cox Stuart test
## 
## data:  Dif_tiempo
## statistic = 11, n = 24, p-value = 0.8388
## alternative hypothesis: non randomness
turning.point.test(Dif_tiempo)
## 
##  Turning Point Test
## 
## data:  Dif_tiempo
## statistic = -0.80561, n = 49, p-value = 0.4205
## alternative hypothesis: non randomness

Los distintos tests de rachas utilizados indican que podemos retener la hipótesis nula, es decir, podemos asumir que los datos son iid.

Procesos de Poisson No Homogéneos (PPNH)

Un Proceso de Poisson No Homogéneo (PPNH) es un proceso en que la tasa de ocurrencia de eventos depende el tiempo. Este tipo de procesos pueden ser, en muchos casos, más realistas que los PP y PPC. Ejemplos de aplicación de PPNH son:

Modelar la llegada de más de cierto número de turistas por mes: claramente habrá meses con alto número de turistas (temporada alta) y meses con menor arribo de turistas (temporada baja).

Modelar el número de retiros de dinero en cajeros automáticos en los días del mes.

¿Otros ejemplos?

Retomando la definición del teórico:

Si N es un proceso de conteo y m es una medida que NO puede expresarse como una constante por la longitud, diremos que N es un Proceso de Poisson No Homogéneo de medida \(m\) (y abreviaremos N es \(PPNH(m)\)) si se cumple:

  1. Para todo intervalo J de los reales positivos, N(J) es una variable aleatoria que tiene distribución de Poisson de parámetro m(J).

  2. Si J, L, M…..es una cantidad arbitraria de intervalos positivos DISJUNTOS, entonces N(J), N(L), N(M),……son variables aleatorias independientes.

Para dejar en claro la diferencia entre los PPNH y los PPC (o el simple PP), recordemos que, en los PPC, los tiempos intereventos son exponenciales de parámetro \(\lambda >0\) e iid. Si N es un PPNH y T1 es el tiempo del primer evento, la distribución de T1 no es exponencial.

Pasemos a simular un proceso de PPNH utilizando la función simpos del paquete {IHSEP}. A partir de esta función lo que haremos es simular seis PP simples (podrían ser más) con distinto lambda y los “uniremos”. Observemos que con esta función tenemos una manera alternativa a la vista más arriba de simular PP simples.

#install.packages("IHSEP")
library(IHSEP)
## Warning: package 'IHSEP' was built under R version 4.3.3
set.seed(69)
tms1 <- simPois(int = function(x) 0.1, cens = 1000) #PP simple con lambda=0.1

T1<-diff(tms1)

#a modo de ejercicio estimo lambda
1/mean(T1) # Lambda = 0.1
## [1] 0.117381
tms2<- simPois(int = function(x) 2, cens = 1000)
#PP simple con lambda=2
T2<-diff(tms2) # Equivalente a T del teórico (Tiempos inter eventos)

Dif_tiempo_NH<- c(T1,T2,T1,T2,T1,T2)

PPNH<- cumsum(Dif_tiempo_NH)

## Visualicemos la cantidad de eventos
hist(PPNH, breaks = 100, main = "Proceso NO homogéneo de Poisson", xlab = "", ylab = "frecuencia")

x<-seq(1, max(PPNH), length=length(PPNH)) 
par(mfrow=c(1,3))
plot(PPNH,x, type="s",ylab="N",xlab="Tiempo")

Vale aclarar que esta función también es útil para simular PPNH directamente si se incluye en la misma una función de intensidad.

Tal como venimos haciendo, asumiendo que no conocemos la estructura de los datos realizaremos test de aleatoriedad y de Lilliefors para determinar qué tipo de proceso siguen nuestros datos.

¿Qué esperamos encontrar?

Aleatoriedad:

#install.packages("randtests")
require(randtests)
runs.test(Dif_tiempo_NH)
## 
##  Runs Test
## 
## data:  Dif_tiempo_NH
## statistic = -3.1642, runs = 3402, n1 = 3534, n2 = 3534, n = 7068,
## p-value = 0.001555
## alternative hypothesis: nonrandomness
cox.stuart.test(Dif_tiempo_NH)
## 
##  Cox Stuart test
## 
## data:  Dif_tiempo_NH
## statistic = 1712, n = 3535, p-value = 0.06428
## alternative hypothesis: non randomness
turning.point.test(Dif_tiempo_NH)
## 
##  Turning Point Test
## 
## data:  Dif_tiempo_NH
## statistic = 1.8994, n = 7071, p-value = 0.05752
## alternative hypothesis: non randomness
#Ajuste a la exponencial de tiempos intereventos
TestLillie=LcKS(Dif_tiempo_NH, "pexp", nreps = 4999)
TestLillie$p.value # devuelve el valor de P
## [1] 2e-04
#install.packages('NHPoisson')
library('NHPoisson')
## Loading required package: stats4
##install.packages("IHSEP")
library(IHSEP)
set.seed(69)
tms <- simPois(int = function(x) 0.1, cens = 1000)

1/mean(diff(tms)) # Lambda = 0.1
## [1] 0.117381
tms2<- simPois(int = function(x) 2, cens = 1000)

PPNH<-c(tms, tms[length(tms)]+tms2, 1996.3+tms, 2992.7+tms2 ) # Poisson no homogeneo

1/mean(diff(PPNH))
## [1] 1.181717
hist(PPNH, breaks = 100, main = "Proceso NO homogéneo de Poisson", 
     xlab = "", ylab = "frecuencia")

Ejercicio

Con la serie de datos enviada, que se basa en los tiempos de los registros de utilizacion de un cajero automático, por favor responda:

  1. Determinar si es un PPH o un PPNH.
  2. Calcular la tasa de incidencia (\(\lambda\)) y el tiempo promedio intereventos (si es PPNH recuerde hacerlo para cada segmento).

Resolvemos 1)

Veamos qué pasa con la independencia (aleatoriedad).

Recordemos lo siguiente:

Comparación entre Proceso de Poisson Homogéneo (PPH) y No Homogéneo (PPNH)
Concepto PPH PPNH
Tasa \(\lambda\) constante \(\lambda(t)\) variable
Intereventos i.i.d. \(\mathrm{Exp}(\lambda)\) Dependientes del tiempo
Incrementos Independientes y estacionarios Independientes, no estacionarios
Diagnóstico Incrementos \(\sim \mathrm{Exp}\) y sin tendencia Tendencia o tasa variable 
path="/Users/lauramontaldo/Desktop/RmdsExtremos/lab2/Datos_LAB2_Poisson.csv"
file.exists(path)
## [1] TRUE
x <- read.csv(path)
head(x)
##            x
## 1 0.04936894
## 2 0.15810106
## 3 0.33718684
## 4 0.38622738
## 5 0.54093769
## 6 0.62420780
# Convertir a vector numérico
x <- x$x

# Verificar
is.numeric(x)
## [1] TRUE
# Gráfico de conteo acumulado N(t)
N <- 0:length(x)
plot(stepfun(x, N),
     main = "Conteo acumulado N(t)",
     xlab = "Tiempo (t)",
     ylab = "Número de eventos")

Es una serie acumulada y estrictamente creciente, porque cada evento ocurre después del anterior.

Definición 3.1 (Proceso de Poisson)

Sea \(\{N(t), \, t \ge 0\}\) un proceso de conteo.
Decimos que es un PPH con tasa \(\lambda > 0\) (denotado \(PP(\lambda)\)) si cumple:

  1. \(N(0) = 0\)
  2. Incrementos independientes: El número de eventos que ocurren en intervalos disjuntos de tiempo son variables aleatorias independientes.
  3. Incrementos estacionarios: El número de eventos en cualquier intervalo de longitud \(t\) tiene distribución de Poisson con media \(\lambda t\), es decir, para todo \(s, t \ge 0\) y \(n = 0, 1, 2, \ldots\) se cumple que: \[ P\big(N(t+s) - N(s) = n\big) = e^{-\lambda t} \frac{(\lambda t)^n}{n!}. \]

Como consecuencia, el proceso cumple que \[ \mathbb{E}[N(t)] = \lambda t, \] por lo que \(\lambda\) se interpreta como la tasa de incidencia o intensidad promedio del proceso.

1- es \(N(0)=0\)?

min(x)
## [1] 0.04936894

Podemos decir que el proceso está cerca de cero.

2- Tiene incrementos independientes: el número de eventos que ocurren en intervalos disjuntos de tiempo son variables aleatorias independientes?

x=sort(x)
incrementos<- diff(x)
plot(incrementos, type = "l")

Al graficar eso da una serie que sube y baja, porque algunos intervalos son más largos y otros más cortos, no tiene por qué ser creciente, ya que cada interevento puede variar aleatoriamente.

  • Prueba de rachas (Runs Test) sobre los incrementos:

Chequamos \(H_0)\) la secuencia es aleatoria vs \(H_1)\) no aleatoriedad en la secuencia. Analizamos si la secuencia de observaciones presenta un orden aleatorio, es decir, si no hay patrón sistemático de aumentos o disminuciones.

 ## Aleatoriedad
library(randtests)
r_test=randtests::runs.test(as.numeric(incrementos))
r_test$p.value
## [1] 0.671844

Como el p-valor es mayor a \(\alpha=0.05\) entonces no \(RH_0\) y concluimos que los datos son aleatorios.

2- Son estacionarios los incrementos? vemos si la secuencia presenta tendencia.

  • Cox–Stuart (tendencia) sobre los incrementos: evalúa si existe una tendencia (creciente o decreciente) en una secuencia de datos en el tiempo.

\(H_0)\) No hay tendencia \(H_1)\) Hay tendencia: los intereventos crecen o decrecen

# Test de tendencia en los tiempos intereventos
cox_res  <- randtests::cox.stuart.test(incrementos, alternative = "two.sided")
print(cox_res)
## 
##  Cox Stuart test
## 
## data:  incrementos
## statistic = 541, n = 1115, p-value = 0.3379
## alternative hypothesis: non randomness

Si no rechazamos \(H_0)\) entonces podemos decir que para un determinado nivel de confianza, hay evidencia de intereventos estacionarios, a tasa constante, y el proceso podría ser PPH. Si rechazamos la \(H_0)\), podemos concluir que hay evidencia, para un determinado nivel de confianza, de que los intereventos crecen y/o decrecen, por lo que la tasa cambia y el proceso podría ser PPNH.

-Turning Point Test (test de puntos de inflexión): otro test de aleatoriedad. Se usa para evaluar si una secuencia de datos (como los tiempos intereventos) muestra dependencia o tendencia local en el tiempo.

Bajo la hipótesis nula \(H_0\):

Si los datos son i.i.d. (sin tendencia ni correlación), el número esperado de turning points es

\[ E[T] = \frac{2}{3}(n - 2) \]

y la varianza está dada por

\[ \mathrm{Var}(T) = \frac{16n - 29}{90}. \]

El estadístico de prueba es aproximadamente normal

\[ Z = \frac{T - E[T]}{\sqrt{\mathrm{Var}(T)}} \sim N(0,1) \]

y el test compara el valor \(Z\) con los límites de una distribución normal estándar.

Testeamos \(H_0)\) La serie es aleatoria (iid) contra \(H_1)\) La serie no es aleatoria.

tp_res <- randtests::turning.point.test(incrementos)
print(tp_res)
## 
##  Turning Point Test
## 
## data:  incrementos
## statistic = 1.7079, n = 2231, p-value = 0.08765
## alternative hypothesis: non randomness

Si el p-valor es mayor a \(\alpha\) entonces no se rechaza \(H_0)\) y hay evidencia de que la serie es aleatoria. Si el p-valor es menor a a \(\alpha\) entonces se rechaza \(H_0)\) y hay evidencia de dependencia o tendencia local.

Los tres tests no muestran evidencia de dependencia o tendencia, para un nivel de confianza del 95\(\%\) por lo que concluimos que son incrementos independientes.

Veamos ahora, si los incrementos \(T_i\), \(i=1,\dots, n\) siguen una distribución exponencial.

Constrastamos \(H_0)\;T_i \sim \text{Exp}(\lambda)\) contra \(H_1)\;T_i\) no es exponencial.

con \(\lambda\) desconocido que se estima como \(\hat{\lambda=1/\bar{T}}\).

require(KScorrect)

TestLillie=LcKS(incrementos, "pexp", nreps = 4999)
TestLillie$p.value # No rechazo H0 es exponencial
## [1] 0.6284

Como no rechazamos \(H_0)\) para un nivel de confianza del 95\(\%\), concluimos que los datos son compatibles con una distribución exponencial, seguimos en la dirección de que el proceso podría llegar a ser PPH.

Veamos ahora graficamente:

qqplot(qexp(ppoints(length(incrementos)), rate = 1/mean(incrementos)),
       incrementos, main = "Q–Q plot: incrementos vs exponencial",
       xlab = "Cuantiles teóricos Exp(hat lambda)", ylab = "Cuantiles observados")
abline(0, 1, col = "red", lwd = 2)

Si los puntos siguen aproximadamente la línea roja, los incrementos son consistentes con una distribución exponencial (iid)

hist(incrementos, breaks = 40, probability = TRUE,
     col = "lightgray", border = "white",
     main = "Distribución de los incrementos",
     xlab = expression(Delta*t))
curve(dexp(x, rate = 1/mean(incrementos)), col = "blue", lwd = 2, add = TRUE)
legend("topright", legend = "Exp(1/mean)", col = "blue", lwd = 2)

Como la curva azul (exponencial ajustada) se superpone bien al histograma, es consistente con un PPH.

Resolvemos 2)

Si la secuencia es

-PPH: una sola tasa y un solo tiempo promedio.

  • PPNH: tasas y tiempos promedio por segmento.

Si el proceso es un , entonces

\[ \widehat{\lambda} = \frac{1}{\overline{T}}, \qquad \overline{T} = \frac{1}{n - 1} \sum_{i = 1}^{n - 1} (x_{i+1} - x_i) \]

# Tasa de incidencia (lambda) y tiempo promedio intereventos
lambda_hat <- 1 / mean(incrementos)
T_promedio <- mean(incrementos)

# Mostrar resultados
cat("Tasa de incidencia (lambda):", lambda_hat, "\n")
## Tasa de incidencia (lambda): 4.465467
cat("Tiempo promedio intereventos:", T_promedio, "\n")
## Tiempo promedio intereventos: 0.2239407