1. Proceso de Poisson homogéneo (HPP)

1.1 Definición

Un proceso de Poisson homogéneo con tasa \(\lambda>0\) es un proceso de conteo \(\{N(t): t\ge 0\}\) tal que:

  1. \(N(0)=0\) y las trayectorias son no-decrecientes por saltos de tamaño 1.
  2. Incrementos independientes: para \(0\le t_0 < t_1 < \cdots < t_k\), los incrementos \(N(t_j)-N(t_{j-1})\) son independientes.
  3. Incrementos estacionarios: \(N(t+s)-N(s)\stackrel{d}{\sim}\text{Poisson}(\lambda t)\).
  4. Pequeños intervalos: para \(h\downarrow 0\), \[ \Pr\{N(t+h)-N(t)=1\}=\lambda h+o(h),\qquad \Pr\{N(t+h)-N(t)\ge 2\}=o(h). \]

Se deduce que: \[ N(t)\sim \text{Poisson}(\lambda t),\qquad \mathbb E[N(t)] = \lambda t,\quad \mathrm{Var}[N(t)]=\lambda t. \]

1.2 Tiempos de arribo e interarribos

Sea \(T_1<T_2<\cdots\) la sucesión de tiempos de arribo y \(X_k=T_k-T_{k-1}\) los interarribos (\(T_0=0\)). Entonces: - \(X_k\stackrel{iid}{\sim}\text{Exp}(\lambda)\). - \(T_n=\sum_{k=1}^n X_k\sim\text{Gamma}(n,\lambda)\) (parámetros: forma \(n\), tasa \(\lambda\)). - Propiedad sin memoria: \(\Pr\{X_k> x+y\mid X_k> x\}=\Pr\{X_k>y\}\).

2. Verosimilitud y MLE (HPP)

2.1 Datos como tiempos de arribo en \([0,T]\)

Observamos \(n\) arribos \(0<T_1<\cdots<T_n\le T\). La verosimilitud (usando independencia de interarribos y c.e.) es \[ L(\lambda\mid T_{1:n},T) =\lambda^n \exp(-\lambda T). \] Equivalente (condicionando en \(N(T)=n\)): los \(T_i\) ordenados son las estadísticas de orden de \(n\) i.i.d. Uniforme\((0,T)\).

El log-verosimil: \[ \ell(\lambda)= n\log\lambda - \lambda T. \]

MLE: \[ \hat\lambda=\frac{n}{T}. \]

Información de Fisher \(I(\lambda)=\frac{T}{\lambda^2}\) \(\Rightarrow\) \(\mathrm{Var}(\hat\lambda)\approx \frac{\lambda^2}{T}\) y por reemplazo \(\mathrm{se}(\hat\lambda)\approx \frac{\hat\lambda}{\sqrt{T}}\).

Un IC (1−α) asintótico para \(\lambda\): \[ \hat\lambda \pm z_{1-\alpha/2}\,\frac{\hat\lambda}{\sqrt{T}}. \]

2.2 Datos como conteos por intervalos

Si solo observas \(Y_j=N(t_j)-N(t_{j-1}) \stackrel{ind}{\sim}\text{Poisson}\big(\lambda \Delta_j\big)\) con \(\Delta_j=t_j-t_{j-1}\), \[ \ell(\lambda)=\sum_{j=1}^m \big[Y_j\log(\lambda \Delta_j)-\lambda \Delta_j - \log(Y_j!)\big] \] \(\Rightarrow\) \(\hat\lambda=\dfrac{\sum_j Y_j}{\sum_j \Delta_j}\).

2.3 Contrastes básicos

  • Wald: \(Z=\dfrac{\hat\lambda-\lambda_0}{\hat\lambda/\sqrt{T}}\approx \mathcal N(0,1)\).
  • Exacto (conteo): si usas \(N(T)\sim\text{Poisson}(\lambda T)\), prueba exacta para \(H_0:\lambda=\lambda_0\).

3. Proceso de Poisson no homogéneo (NHPP)

Ahora la tasa es función del tiempo \(\lambda(t)>0\). Defina el compensador \[ \Lambda(t)=\int_0^t \lambda(s)\,ds. \] Entonces: \[ N(t)\sim \text{Poisson}\big(\Lambda(t)\big),\qquad \Pr\{N(t+h)-N(t)=1\}=\lambda(t)\,h+o(h). \]

Si observamos tiempos de arribo \(0<t_1<\cdots<t_n\le T\), la verosimilitud es \[ L(\theta)\;=\;\prod_{i=1}^n \lambda(t_i;\theta)\;\exp\!\big(-\Lambda(T;\theta)\big),\quad \ell(\theta)=\sum_{i=1}^n \log \lambda(t_i;\theta)-\Lambda(T;\theta). \] Ejemplo paramétrico: \(\lambda(t;\theta)=\exp(\theta_0+\theta_1 t)\) (siempre positiva). Entonces \(\Lambda(T;\theta)=\int_0^T \exp(\theta_0+\theta_1 s)ds =\frac{e^{\theta_0}}{\theta_1}\big(e^{\theta_1 T}-1\big)\) (si \(\theta_1\neq 0\)); si \(\theta_1=0\), \(\Lambda(T)=e^{\theta_0}T\).

4. Simulación

4.1 HPP por interarribos

set.seed(123)
lambda <- 2      # eventos por unidad de tiempo
Tend   <- 10     # horizonte

# Simular interarribos ~ Exp(lambda) hasta exceder T
x <- rexp(10000, rate = lambda)
t <- cumsum(x)
t <- t[t <= Tend]
n <- length(t)
n
## [1] 23
# Conteo total en [0, T]
N_T <- rpois(1, lambda = lambda * Tend)
N_T
## [1] 27
set.seed(321)
lambda_fun <- function(tt) exp(0.1 + 0.2*tt)          # ejemplo NHPP
Tend <- 5
lam_max <- lambda_fun(Tend)                           # cota sencilla

# Simulación por thinning
t <- c(); s <- 0
while(TRUE){
  s <- s + rexp(1, rate = lam_max)                   # arribo candidato
  if(s > Tend) break
  if(runif(1) < lambda_fun(s)/lam_max) t <- c(t, s)  # aceptar
}
length(t); head(t)
## [1] 12
## [1] 0.05496428 0.47278242 0.84078986 1.01810261 1.07946531 1.16783462
# Dados tiempos t en [0,T], el MLE es n/T
Tend <- 10
lambda_hat <- length(t)/Tend
se_hat     <- lambda_hat / sqrt(Tend)
c(lambda_hat = lambda_hat, se = se_hat,
  CI95_L = lambda_hat - 1.96*se_hat,
  CI95_U = lambda_hat + 1.96*se_hat)
## lambda_hat         se     CI95_L     CI95_U 
##  1.2000000  0.3794733  0.4562323  1.9437677
set.seed(42)
breaks <- c(0, 1, 1.5, 3, 5, 6.2, 10)
Delta  <- diff(breaks)
Y      <- rpois(length(Delta), lambda = 1.8 * Delta)   # simulado

lambda_hat <- sum(Y) / sum(Delta)
lambda_hat
## [1] 2.3
# Suponiendo que tenemos tiempos de arribo 't' en [0, Tend]
lambda_log <- function(tt, theta) exp(theta[1] + theta[2]*tt)
Lambda_log <- function(Tend, theta){
  if(abs(theta[2]) < 1e-10) return( exp(theta[1]) * Tend )
  exp(theta[1]) * (exp(theta[2]*Tend) - 1) / theta[2]
}

loglik <- function(theta, t, Tend){
  sum(log(lambda_log(t, theta))) - Lambda_log(Tend, theta)
}
theta0 <- c(log(length(t)/Tend), 0)  # inicial: tasa const.
fit <- optim(theta0, fn = function(th) -loglik(th, t=t, Tend=Tend),
             control = list(reltol = 1e-10), hessian = TRUE, method = "BFGS")
fit$par
## [1]  2.0207267 -0.6274665
# Varianza asintótica (inversa de -Hessiano)
vcov <- solve(fit$hessian)
se   <- sqrt(diag(vcov))
cbind(estimate = fit$par, se = se)
##        estimate        se
## [1,]  2.0207267 0.4138214
## [2,] -0.6274665 0.1882754
# Diagnóstico HPP
lambda_hat <- length(t)/Tend
x_hat <- diff(c(0, t))                   # interarribos
# QQ-plot exponencial
ExpQuant <- function(p, rate) -log(1-p)/rate
pp <- (1:length(x_hat) - 0.5)/length(x_hat)
plot(sort(x_hat), ExpQuant(pp, rate = lambda_hat),
     xlab="Interarribos observados (ordenados)",
     ylab="Cuantiles Exp(lambda_hat)",
     main="QQ-plot Exponencial vs interarribos")
abline(0, 1)

set.seed(7)
lambda_true <- 3; Tend <- 8
x <- rexp(10000, rate=lambda_true)
t_e1 <- cumsum(x); t_e1 <- t_e1[t_e1 <= Tend]
lambda_hat <- length(t_e1)/Tend
se_hat     <- lambda_hat/sqrt(Tend)
IC <- lambda_hat + c(-1,1)*1.96*se_hat
c(lambda_hat=lambda_hat, IC_L=IC[1], IC_U=IC[2])
## lambda_hat       IC_L       IC_U 
##  2.5000000  0.7675884  4.2324116
Tend <- 10; Nobs <- 14; lambda0 <- 1.2
mu0 <- lambda0 * Tend
pval <- 2*min(ppois(Nobs, mu0), 1-ppois(Nobs-1, mu0))
c(p_value = pval)
##   p_value 
## 0.6369287
set.seed(99)
lambda_fun <- function(tt) exp(0.2 + 0.3*tt)
Tend <- 5; lam_max <- lambda_fun(Tend)
t <- c(); s <- 0
while(TRUE){
  s <- s + rexp(1, rate = lam_max)
  if(s > Tend) break
  if(runif(1) < lambda_fun(s)/lam_max) t <- c(t, s)
}

lambda_log <- function(tt, theta) exp(theta[1] + theta[2]*tt)
Lambda_log <- function(Tend, theta){
  if(abs(theta[2]) < 1e-10) return( exp(theta[1]) * Tend )
  exp(theta[1]) * (exp(theta[2]*Tend) - 1) / theta[2]
}
loglik <- function(theta, t, Tend){
  sum(log(lambda_log(t, theta))) - Lambda_log(Tend, theta)
}
theta0 <- c(log(length(t)/Tend), 0)
fit <- optim(theta0, fn = function(th) -loglik(th, t=t, Tend=Tend),
             method="BFGS", hessian=TRUE)
vcov <- solve(fit$hessian); se <- sqrt(diag(vcov))
rbind(estimate = fit$par, se = se)
##               [,1]      [,2]
## estimate 0.2820649 0.1879404
## se       0.6859498 0.2134965
theta_hat <- fit$par
u <- sapply(t, function(tt) Lambda_log(tt, theta_hat))
du <- diff(c(0, u))
pp <- (1:length(du)-0.5)/length(du)
plot(sort(du), -log(1-pp), xlab="du (ordenados)", ylab="Cuantiles Exp(1)",
     main="QQ-plot Exp(1) de incrementos transformados")
abline(0,1)