Modelo de Poisson con cambio de régimen
Alfredo Zea Valenzuela
09/12/2022
1. Modelo de Poisson con cambio de régimen
Un Modelo de Poisson con cambio de régimen se puede expresar de la siguiente manera: \[ y_{i} \sim Poisson(\lambda_{1})~~~~si~~i=1,2,...,m \] \[ y_{i} \sim Poisson(\lambda_{2})~~~~si~~i=m +1,..,N \] Se advierte que las primeras m observaciones tienen una distribución \(Poisson\) con parámetro \(\lambda_{1}\) y a partir de la observación (m+1) sigue una distribución de \(Poisson\) con parámetro \(\lambda_{2}\), en resumen, a partir de la observación (m+1) se produce un cambio de régimen.
2. Simulación de una muestra
A continuación se procede a simular una muestra con cambió del régimen con parámetros: \(\lambda_{1}\),\(\lambda_{2}\) y \(m\) :
# Eliminar los objetos de Global Environment
rm(list = ls())
options("digits" = 4)
set.seed(987)
# ==============================================================================
# Simulando un modelo de Poisson con cambio de régimen
# ==============================================================================
# Total de periodos
N = 50
# Periodo de cambio de regimen
m = 20
# Primer regimen
# Yi ~ Poison(5) i: 1..m
λ1 = 3
Y1 = rpois(m,λ1)
media1 = round(mean(Y1),3)
# Segundo regimen
# Yi ~ Poison(15) i: 1..m
λ2 = 7
Y2 = rpois( (N-m ),λ2)
media2 = round(mean(Y2),3)
df_resumen =
data.frame( Parametro = c('λ1','λ2','m'),
Teorico = c(λ1,λ2,m),
MediaMV = c(media1,media2,''))
df_resumen## Parametro Teorico MediaMV
## 1 λ1 3 3
## 2 λ2 7 6.8
## 3 m 20# Regimen 1 unión Regimen 2
Y = c(Y1,Y2)
par(mfrow=c(2,2))
hist(Y1,main = sprintf("Regimen 1: Poison(%.0f)",λ1),
ylab = "Frecuencia",
xlab = 'y', breaks = 20)
abline(v=media1, col='green', lwd=3)
hist(Y2,main = sprintf("Regimen 2: Poison(%.0f)",λ2),
ylab = "Frecuencia",
xlab = 'y', breaks = 20)
abline(v=media2, col='red', lwd=3)
hist(Y,main = "Regimen 1 y Regimen 2",
ylab = "Frecuencia",
xlab = 'y', breaks = 20)
abline(v=media1, col='green', lwd=3)
abline(v=media2, col='red', lwd=3)
par(mfrow=c(1,1))
Mediante el enfoque clásico no es posible estimar el parámetro m que corresponde a la posición de la observación donde hay un cambio de régimen.
3. Modelo con cambio de régimen con enfoque bayesiano
Aplicando el enfoque bayesiano, extendemos el modelo inicial considerando que \(\lambda_{1}\), \(\lambda_{2}\) y \(m\) son variables aleatorias y por tanto tienen una distribución a priori:
\[ y_{i} \sim Poisson(\lambda_{1})~~con~\lambda_{1}>0~~~~si~~i=1,2,...,m \] \[ y_{i} \sim Poisson(\lambda_{2})~~con~\lambda_{2}>0~~~~si~~i=m +1,..,N \] \[ \lambda_{1} \sim Gamma(a,b)~~con~a>0~y~b>0\] \[ \lambda_{2} \sim Gamma(c,d)~~con~c>0~y~d>0\] \[ m \sim Uniforme\{1,2,...N\}\]
En este modelo tambien se incluyen los Hiperparámetros: a,b,c y d.
4. Distribuciones condicionales de los parámetros
Si deseamos aplicar el algoritmo de Gibbs necesitamos contar con las distribuciones condicionales de cada uno de los parámetros.
Para el parámetro \(\lambda_{1}\) se puede determinar analíticamente la distribución condicional dado que se utilizó un priori conjugado:
\[\lambda_{1}/\lambda_{2},m \sim Gamma(\sum_{i=1}^{m}y_{i} + a , b + m)\] De manera similar, la distribución condicional de \(\lambda_{2}\) es:
\[\lambda_{2}/\lambda_{1},m \sim Gamma(\sum_{i=m+1}^{N}y_{i} + c , N - m +d)\] Para el caso del parámetro \(m\), no es posible determinar analíticamente la distrución condicional. Sin embargo, se cuenta con el kernell de la distribución condicional:
\[\pi(m/\lambda_{1},\lambda_{2}) \propto \lambda_{1}^{\sum_{i=1}^{m}y_{i}}e^{-m\lambda_{1}}\lambda_{2}^{\sum_{i=m+1}^{N}y_{i}}e^{m\lambda_{2}}\] A continuación se muestra un conjunto funciones que implementan la obtención aleatoria de un valor de \(m\) dado los valores de \(\lambda_{1}\) y \(\lambda_{2}\) utilizando el método de simulación de tranformación inversa para obtener una muestra aleatoria para el caso de distribuciones discretas:
fKProbabilidad: Evalua el Kernell de la distribución \(\pi\) dados de \(\lambda_{1}\), \(\lambda_{2}\) y \(m\):
# ------------------------------------------------------------------------------
# Función del kernel de probabilidad de m/λ1,λ2
# fKProbabilidad(Y,2,5,15)
# ------------------------------------------------------------------------------
fKProbabilidad = function(pY,plambda1,plambda2,pm){
lN = length(pY)
lsuma1 = sum(pY[1:pm])
if (pm < lN){
lsuma2 = sum(pY[(pm +1):lN])
}
else {
lsuma2 = 0.00
}
llt1 = lsuma1*log(plambda1)
llt2 = -pm*plambda1
llt3 = lsuma2*log(plambda2)
llt4 = pm*plambda2
lkp = llt1 + llt2 + llt3 + llt4
return( exp(lkp))
}
fKProbabilidad(Y,2,5,15)## [1] 3.647e+186fDistribucionProbabilidad: Calcula para cada \(m\) dados \(\lambda_{1}\) y \(\lambda_{2}\) un valor proporcional a su probabilidad \(kP\):
# ------------------------------------------------------------------------------
# Función del kernel de distribución de probabilidad de m/λ1,λ2
# fDistribucionProbabilidad(Y,12,5)
# ------------------------------------------------------------------------------
fDistribucionProbabilidad = function( pY,plambda1,plambda2){
lN = length(pY)
MDistribucion = matrix( 0.00 ,lN,2)
for (i in 1:lN) {
MDistribucion[i,1] = i
MDistribucion[i,2] = fKProbabilidad(pY,plambda1,plambda2,i)
}
colnames(MDistribucion) = c("m","kP")
return(MDistribucion)
}
fDistribucionProbabilidad(Y,12,5)## m kP
## [1,] 1 4.253e+182
## [2,] 2 4.269e+182
## [3,] 3 5.381e+180
## [4,] 4 6.783e+178
## [5,] 5 2.052e+177
## [6,] 6 1.078e+175
## [7,] 7 1.359e+173
## [8,] 8 7.137e+170
## [9,] 9 5.182e+169
## [10,] 10 1.568e+168
## [11,] 11 3.431e+165
## [12,] 12 1.802e+163
## [13,] 13 3.944e+160
## [14,] 14 2.072e+158
## [15,] 15 4.534e+155
## [16,] 16 5.715e+153
## [17,] 17 1.251e+151
## [18,] 18 5.231e+150
## [19,] 19 1.145e+148
## [20,] 20 3.464e+146
## [21,] 21 8.344e+146
## [22,] 22 3.490e+146
## [23,] 23 2.534e+145
## [24,] 24 6.104e+145
## [25,] 25 2.033e+147
## [26,] 26 1.068e+145
## [27,] 27 4.466e+144
## [28,] 28 1.868e+144
## [29,] 29 3.255e+143
## [30,] 30 2.363e+142
## [31,] 31 9.883e+141
## [32,] 32 9.921e+141
## [33,] 33 1.729e+141
## [34,] 34 1.735e+141
## [35,] 35 5.250e+139
## [36,] 36 3.812e+138
## [37,] 37 6.643e+137
## [38,] 38 2.778e+137
## [39,] 39 1.606e+138
## [40,] 40 9.287e+138
## [41,] 41 3.093e+140
## [42,] 42 3.104e+140
## [43,] 43 2.254e+139
## [44,] 44 1.637e+138
## [45,] 45 2.063e+136
## [46,] 46 1.498e+135
## [47,] 47 3.609e+135
## [48,] 48 4.549e+133
## [49,] 49 1.903e+133
## [50,] 50 7.958e+132frDistribucionProbabilidad : Genera un valor aleatorio de \(m\) dados \(\lambda_{1}\) y \(\lambda_{2}\) :
# ------------------------------------------------------------------------------
# Función que genera un m aletorio de la distribución probabilidad de m/λ1,λ2
# frDistribucionProbabilidad(Y,12,5)
# ------------------------------------------------------------------------------
frDistribucionProbabilidad = function( pY,plambda1,plambda2){
D = fDistribucionProbabilidad(pY,plambda1,plambda2)
lm = sample(x = D[,"m"], 1, replace = T, prob =D[,"kP"])
return( lm )
}
frDistribucionProbabilidad(Y,12,5)## [1] 2Con estas funciones de disribuciones condicionales para generar valores aleatorios de \(m\), \(\lambda_{1}\) y \(\lambda_{2}\) es posible implementar el algoritmo de muestreo de Gibbs.
5. Aplicación del algoritmo de Gibbs
A continuación se muestra la función de implementa el algoritmo de Gibbs para obtener una muestra de una distribución conjunta \(\pi(\lambda_{1},\lambda_{2},m | Y)\) conociendo sus distribuciones condicionales.
# ------------------------------------------------------------------------------
# Función de implementa el algoritmo de Gibs
# ------------------------------------------------------------------------------
fMuestradorGibbs = function( pY,pP,pQ,pOi,pHi){
a = pHi[1]
b = pHi[2]
c = pHi[3]
d = pHi[4]
Os = matrix(0.00, nrow = pP , ncol = length(pOi))
Os[1,] = pOi
N = length(pY)
for (i in 2:pP){
suma1 = sum(pY[1:Os[(i-1),3 ]])
if (Os[(i-1),3] < N) {
suma2 = sum(pY[(Os[(i-1),3] +1):N])
} else {
suma2 = 0.00
}
Os[i,1] = rgamma(n = 1 ,suma1 + a , b + Os[(i-1),3] )
Os[i,2] = rgamma(n = 1 ,suma2 + c , N - Os[(i-1),3] + d)
Os[i,3] = frDistribucionProbabilidad( pY,Os[i,1],Os[i,2])
if (i%%floor(pP/10) == 0){
cat(100*round(i/pP, 1), "% completado ... \n", sep = "" )
}
}
return(Os[(pP - pQ + 1):pP,])
}Procedemos a simular la muestra de los parámetros \(\lambda_{1}\) y \(\lambda_{2}\) y \(m\) con el objetivo de calcular media, mediana y moda de sus distribuciones.
# ==============================================================================
# Obtención de muestra de párametros mediante del algoritmo de Gibbs
# ==============================================================================
# ------------------------------------------------------------------------------
# Simulación
# ------------------------------------------------------------------------------
# Número de iteraciones
P= 100000
# Tamaño de la muestra
Q = 10000
# Vector de parametros inicial λ1,λ2,m
Oi = c(1,1,1)
# Vector de Hiperparámetros
Hi = c(0.1,0.1,0.1,0.1)
Os = fMuestradorGibbs(Y,P,Q,Oi,Hi)## 10% completado ...
## 20% completado ...
## 30% completado ...
## 40% completado ...
## 50% completado ...
## 60% completado ...
## 70% completado ...
## 80% completado ...
## 90% completado ...
## 100% completado ...tail(Os,50)## [,1] [,2] [,3]
## [9951,] 3.477 6.264 17
## [9952,] 3.087 6.551 19
## [9953,] 2.211 7.088 19
## [9954,] 2.624 6.519 20
## [9955,] 2.849 6.736 19
## [9956,] 3.050 6.765 19
## [9957,] 3.481 7.219 17
## [9958,] 2.439 6.190 17
## [9959,] 2.690 6.157 20
## [9960,] 2.805 6.774 20
## [9961,] 2.911 6.962 20
## [9962,] 2.759 6.754 17
## [9963,] 2.759 6.244 20
## [9964,] 3.480 5.881 20
## [9965,] 3.570 6.971 20
## [9966,] 2.789 5.757 16
## [9967,] 2.892 6.620 17
## [9968,] 3.021 6.623 20
## [9969,] 2.794 7.646 20
## [9970,] 2.760 7.529 20
## [9971,] 2.989 6.126 18
## [9972,] 3.317 6.031 19
## [9973,] 2.698 6.944 20
## [9974,] 3.456 7.018 19
## [9975,] 2.888 6.449 20
## [9976,] 3.056 6.361 20
## [9977,] 2.447 6.990 20
## [9978,] 3.485 6.441 17
## [9979,] 3.279 6.836 20
## [9980,] 2.919 6.930 20
## [9981,] 2.673 6.975 20
## [9982,] 3.029 6.508 17
## [9983,] 3.047 6.686 20
## [9984,] 3.604 6.805 17
## [9985,] 2.153 6.640 19
## [9986,] 2.798 6.835 16
## [9987,] 2.498 6.148 17
## [9988,] 2.530 5.970 19
## [9989,] 3.037 6.461 21
## [9990,] 3.741 6.238 20
## [9991,] 2.698 6.908 20
## [9992,] 3.372 6.654 19
## [9993,] 3.870 6.620 20
## [9994,] 3.029 7.130 20
## [9995,] 2.743 6.288 20
## [9996,] 2.831 7.282 19
## [9997,] 2.986 6.991 19
## [9998,] 2.833 6.139 19
## [9999,] 2.971 6.186 20
## [10000,] 2.862 6.609 20Se muestra los ultimos 50 valores de la muestra.
6. Estimación puntual de los parámetros
# ------------------------------------------------------------------------------
# Función que obtienen la moda de un vector de datos
# ------------------------------------------------------------------------------
fObtenerModa <- function(v) {
uniqv <- unique(v)
uniqv[which.max(tabulate(match(v, uniqv)))]
}
# ------------------------------------------------------------------------------
# Analisis descriptivo de parámetro
# ------------------------------------------------------------------------------
fAnalisisDescriptivo = function( Os,s){
ls_title = "Distribucion de parámetro "
xmedia = mean (Os)
xmedian = median(Os)
xmoda = fObtenerModa(Os)
xvar = var(Os)
cat("Estadisticos de la muestra \n",
"media :" ,xmedia,"\n",
"mediana :" ,xmedian,"\n",
"moda :" ,xmoda,"\n",
"varianza:" ,xvar)
par(mfrow=c(2,2))
plot( Os,type="l",xlab = "Iteraciones" , ylab = expression(beta),
main = "Convergencia de parámetro")
acf(Os,main = "Autocorrelación de la muestra")
hist(Os,main = ls_title,
ylab = "Frecuencia",
xlab = s, prob = TRUE, breaks = 50)
abline(v=xmedia, col='blue', lwd=3, lty=2)
par(mfrow=c(1,1))
return(c(xmedia,xmedian,xmoda,xvar))
}Parámetro: \(\lambda_{1}\)
parm1= fAnalisisDescriptivo(Os[,1],"λ1")## Estadisticos de la muestra
## media : 2.961
## mediana : 2.945
## moda : 3.524
## varianza: 0.1627
Parámetro: \(\lambda_{2}\)
parm2= fAnalisisDescriptivo(Os[,2],"λ2")## Estadisticos de la muestra
## media : 6.704
## mediana : 6.69
## moda : 7.033
## varianza: 0.2299
Parámetro: \(m\)
parm3= fAnalisisDescriptivo(Os[,3],"m")## Estadisticos de la muestra
## media : 19.21
## mediana : 20
## moda : 20
## varianza: 1.363
7. Cuadro resumen de la estimación de los parámetros
A continuación se muestra un cuadro resumen de los resultados obtenidos de la estimación:
df_resumen[,"MediaBY"] = c(parm1[1],parm2[1],parm3[1])
df_resumen[,"MedianaBY"] = c(parm1[2],parm2[2],parm3[2])
df_resumen[,"ModaBY"] = c(parm1[3],parm2[3],parm3[3])
df_resumen[,"VarBY"] = c(parm1[4],parm2[4],parm3[4])
df_resumen## Parametro Teorico MediaMV MediaBY MedianaBY ModaBY VarBY
## 1 λ1 3 3 2.961 2.945 3.524 0.1627
## 2 λ2 7 6.8 6.704 6.690 7.033 0.2299
## 3 m 20 19.215 20.000 20.000 1.3626