Introducción.

En esta publicación se muestra como se estiman los parámetros de un modelo de Binomial- Beta con punto de cambio.

Modelo Binomial-Beta con punto de cambio.

Suponga que en un casino a diario se rifan 10 premios de 1.000 dolares entre los clientes. El gerente del casino asegura que la probabilidad de que un premio quede en poder del público es del \(30\%\). Peter además de ser un apostador empedernido, es amante de la estadística y considera que la probabilidad de ganar uno de los 10 premios rifados a diario a caído al \(10\%\). Peter decide estimar el punto de cambio y la probabilidad de ganar antes y después de la caida en la probabilidad. Para lograr esto cuenta con los siguientes datos:

set.seed(1036944)
y <- c(rbinom(n = 30, size = 10, prob=0.3),rbinom(n=30, size=10 , prob=.1))
y
##  [1] 1 4 1 0 6 2 5 6 3 3 1 4 1 3 3 3 4 1 5 6 1 5 2 4 2 5 2 2 2 4 1 0 1 0 0
## [36] 2 1 1 0 2 2 2 1 1 0 2 0 0 0 1 1 0 0 3 2 0 1 2 0 0

Como se puede observar los primeros 30 datos son generados desde una distribución binomial con parámetros \(\psi=0.3\) y \(n=10\) y los últimos 30 se generan desde una distribución binomial con \(\theta=0.1\) y \(n=10\). Asumiendo k desconocido, simbólicamente se tiene que: \[y_i \sim binom(n=10,\psi=0.3)\ \ \ \text{para} \ \ i=\,1,\,2,\,\ldots\,k \] \[y_i \sim binom(n=10,\theta=0.1)\ \ \ \text{para} \ \ i=\,k+1,\,k+2,\,\ldots\,N \]

La función de verosimilitud en este caso está dada por:

\[\begin{align} L\left(\psi,\,\theta,\,k\,|Datos\right)&\propto\prod^k_{i=1}\psi^{y_i}(1-\psi)^{n-y_i}\prod^{N}_{i=k+1}\theta^{y_i}(1-\theta)^{n-y_i}\\ &\propto \psi^{\sum^{k}_{i=1}y_i}(1-\psi)^{nk-\sum^{k}_{i=1}y_i}\theta^{\sum^{N}_{i=k+1}y_i}(1-\theta)^{(N-k)n-\sum^{N}_{i=k+1}y_i} \end{align}\]

Las funciones de distribución a priori están dadas por: \[p(k)=1/N \, \, \, \text{para} \, \, \, k=\,1,\,2\,\ldots\,N\] \[p(\psi)\propto \psi^{a_1-1}(1-\psi)^{b_1-1}\] \[p(\theta)\propto \theta^{a_2-1}(1-\theta)^{b_2-1}\]

La distribución conjunta a priori para \(\theta\), \(\psi\), y \(k\) está dada por \(p(\psi,\,\theta,\,k)=p(\psi)\,p(\theta)\,p(k)\). La distribución a posteriori está dada por: $$ \[\begin{align} p(\psi,\,\theta,\,k|Datos)&\propto L\left(\psi,\,\theta,\,k\,|Datos\right)p(\psi,\,\theta,\,k)\\ &\propto \psi^{\sum^{k}_{i=1}y_i+a_1-1}(1-\psi)^{nk-\sum^{k}_{i=1}y_i+b_1-1}\theta^{\sum^{N}_{i=k+1}y_i+a_2-1}(1-\theta)^{(N-k)n-\sum^{N}_{i=k+1}y_i+b_2-1} \end{align}\] $$ Donde \(a_1 = 0.4\), \(b_1 = 3.6\), \(a_2 = 1.2\) y \(b_2 = 2.8\). Para calcular el estimador bayesiano dado por la media de los parámetros de la distribución a posterior se utiliza el muestreador de gibbs. Por tanto se definen las distribuciones condicionales a posteriri para \(\psi\), \(\theta\) y \(k\):

\[p(\psi| \, \theta,\,k,\,Datos)\propto \psi^{\sum^{k}_{i=1}y_i+a_1-1}(1-\psi)^{nk-\sum^{k}_{i=1}y_i+b_1-1} \] \[p(\theta| \, \psi,\,k,\,Datos)\propto\theta^{\sum^{N}_{i=k+1}y_i+a_2-1}(1-\theta)^{(N-k)n-\sum^{N}_{i=k+1}y_i+b_2-1} \] \[p(k| \, \psi,\,\theta,\,Datos)\propto \left(\frac{\psi}{\theta}\right)^{\sum^{k}_{i=1}y_i}\left(\frac{1-\psi}{1-\theta}\right)^{nk-\sum^{k}_{i=1}y_i} \]

Aplicación del Muestrador de Gibbs y estimación

El muestrador de Gibbs se aplica de la siguiente manera:

n <- length(y) #length of the data
m <- 10000 #length of the chain
p1 <- p2 <- k <- numeric(m)
L <- numeric(n)
k[1] <- sample(1:n, 1)
p1[1] <- 0.5
p2[1] <- 0.5
a1 <- 0.4
b1 <- 3.6
a2 <- 1.2
b2 <- 2.8 
# run the Gibbs sampler
for (i in 2:m) {
  kt <- k[i-1]
  #generate mu
  r1 <- sum(y[1:kt])+a1
  r2 <-  kt*10-sum(y[1:kt])+b1
  p1[i] <- rbeta(1, shape1 = r1, shape2 = r2)
  #generate lambda
  if (kt + 1 > n){
    r3 <- a2 
    r4 <- b2}else{
      r3 <- sum(y[(kt+1):n])+a2
      r4 <- (n-kt)*10-sum(y[(kt+1):n])+b2 
    }
  p2[i] <- rbeta(1, shape1 = r3, shape2 = r4)
  
  for (j in 1:n) {
    L[j] <- ((1-p1[i])/(1-p2[i]))^(10*j-sum(y[1:j]))*(p1[i]/p2[i])^(sum(y[1:j]))  }
  L <- L / sum(L)
  #generate k from discrete distribution L on 1:n
  k[i] <- sample(1:n, prob=L, size=1)
}

Los estimadores bayesianos son el promedio de cada una de las muestras generadas para cada parámetro y están dados por:

mean(k[-c(1:1000)])
## [1] 30.05656
mean(p1[-c(1:1000)])
## [1] 0.2992543
mean(p2[-c(1:1000)])
## [1] 0.0906856

Para \(k\), \(p_1\) y \(p_2\) respectivamente. Se puede observar que los estimadores obtenidos son cercanos a los verdaderos valores de los parámetros