Estimadores bayesianos de una proporción

Tenemos que $X \sim Bern(p)$ y la distribución a priori de p es $p \sim Beta(a,b)$. Por tanto, la distribución a posteriori es, $p|\textbf{x} \sim Beta( a+n\widehat{p} ,b+n-n\widehat{p})$

Veamos esto último. La función de probabilidad marginal de p es $\xi(p) = \frac{1}{B(a,b)} \cdot p^{a-1}(1-p^{b-1})$. La función de probabilidad de X es $ f_p(x) = p^x(1-p){1-x} $ y la función de verosimilitud será \[ L(p|\textbf{x}) = f_n(\textbf{x}|p) = \prod_{i=1}^n f_p(x) = \prod_{i=1}^n p^{x_i}(1-p)^{1-x_i} = p^{\sum_{i=1}^\infty x_i} (1-p)^{n-\sum_{i=1}^\infty x_i} = p^{n\widehat{p}} (1-p)^{n-n\widehat{p}} = p^k (1-p)^{n-k}\] Por tanto, como el logartimo es una función monótona creciente equivale a maximizar $\lambda(p|\textbf{x}) = \log(L(p|\textbf{x})) = k\log(p) + (n-k)\log(1-p)$ y derivando nos queda $\frac{\partial{\lambda(p|\textbf{x})}}{\partial{p}} = \frac{k}{p} + \frac{n-k}{p} = 0 \iff p=\frac{k}{n}$ luego el MLE es la proporción muestral.

Usando la fórmula de de estimación de Bayes la función de probabilidad a posteriori es \[\xi(p|\textbf{x}) = \frac{\xi(p)f_n(\textbf{x}|p)}{f_n(\textbf{x})} = \frac{1}{B(a,b)f_n(\textbf{x})}p^{a+k-1}(1-p)^{n-k+b-1}\] Por tanto, $p|\textbf{x} \sim Beta(a+k,b+n-k)$ es la distribución a posteriori. Por tanto, el estimador bayesiano es \[p_b = \frac{a+k}{a+b+n}\]

Tenemos que las ventajas del estimador bayesiano son sus dotes de estimar el parámetro según su distribución de probabilidad a posteriori.

Ejercicio 1

Veamos experimentalmente como el estimador bayesiano nos da una muy buena aproximación de la esperanza de la distribución a posteriori. De hecho, es mucho más fina que la del MLE.

ebc<-function(a,b,k,n,L){
pb<-(a+k)/(a+b+n)    # estimador bayesiano
vp<-c()
for(i in 1:L){
  u<-0
  while(u!=k){
    p<-rbeta(1,a,b)  # tomamos el parámetro a priori
    u<-rbinom(1,n,p) # hallamos el valor para la variable supuesto el parámetro
  }                  # cuando u es el número de aciertos que ha habido en la muestra estamos en el caso deseado
  vp[i]<-p           # entonces guardamos p como la probabilidad a posteriori que nos ha dado la muestra que queríamos
}
print(pb)            # Imprimimos el estimador bayesiano
print(k/n)           # Imprimimos el MLE
return(mean(vp))     # Imprimimos el promedio de los valores deseados
}
ebc(3,4,4,13,10000)

## [1] 0.35
## [1] 0.3076923

## [1] 0.3493133

Ejercicio 2

Comprobamos ahora la propiedad de que el estimador bayesiano minimiza el error cuadrático medio. Lo vemos dibujando una gráfica del error cuadrático medio y comprobando que en el mínimo se alcanza ese valor. También vemos que el valor de la varianza de la distribución p a posteriori coincide con ese valor.

ecmg<-function(a,b,k,n,L){
  pb<-(a+k)/(a+b+n)   # estimador bayesiano
  pm<-k/n             # MLE
  vp<-c()             # idem construcción
  for(i in 1:L){
    u<-0
    while(u!=k){
      p<-rbeta(1,a,b)
      u<-rbinom(1,n,p)
    }
    vp[i]<-p
  }
  puntos<-seq(0,1,0.01)                                 # una partición del intervalo (0,1) con norma 0.01
  ecm<-function(t)mean((vp-t)^2)                        # función error cuadrático medio
  valores<-sapply(puntos,ecm)                           # aplicamos la función a los puntos de la partición
  plot(puntos,valores,type='l',xlab='x',ylab='ecm')     # esbozamos una curva con puntos en abscisas y valores en ordenadas
  segments(pm,0,pm,0.1,col='red')                       # Abscisa del MLE
  segments(pb,0,pb,0.1,col='green')                     # Abscisa del bayesiano
  print(ecm(pb))                                        # error cuadrático medio en el punto pb
  print((a+k)*(b+n-k)/(a+b+n)^2/(a+b+n+1))              # la varianza de la distribución a posteriori (beta)
  print(var(vp))                                        # varianza de los valores deseados
}
ecmg(3,4,4,13,10000)

## [1] 0.01089505
## [1] 0.01083333
## [1] 0.01089508

Ejercicio 3

Veamos cuando la proporción muestral es igual a la esperanza de la distribución a priori.

Tenemos que la distribución a priori de p es $p \sim Beta(a,b)$, luego su esperanza a priori es $\xi(p) = \frac{a}{a+b}$. Sabemos que tiene que ser igual a la proporción muestral, esto es, $\frac{a}{a+b} = \frac{k}{n} = \widehat{p}$. Despejando $k$ tenemos que $k = \frac{an}{a+b}$ y sustituyendo esto en la esperanza de la distribución a posteriori nos queda \[E(p|\textbf{x}) = \frac{a+k}{a+n+b} = \frac{a+\frac{an}{a+b}}{a+n+b} = \frac{a^2+ba+an}{(a+b)(a+n+b)} = \frac{a(a+n+b)}{(a+b)(a+n+b)} = \frac{a}{a+b}\] luego el estimador bayesiano coincide con el estimador de máxima verosimilitud.

Veamos que el bayesiano solo coinde con el MLE en este caso. Partimos de que la esperanza de la distribución a posteriori,i.e. estimador bayesiano, es $E(p|\textbf{x}) = \frac{a+k}{a+b+n}$ y la esperanza de la distribución a priori es $E(p) = \frac{a}{a+b}$, e igualando estas expresiones nos queda \[ \frac{a+k}{a+b+n}=\frac{a}{a+b} \iff a^2+ab+ka+kb = a^2+ab+an \iff k(a+b)=an \iff \frac{a}{a+b} = \frac{k}{n} \] luego solo sucede en el caso en el que la esperanza a priori coincida con la proporción muestral.

En una distribución geométrica

Ejercicio 4

La distribución a priori de p es $p \sim Unif(0,1)$ luego $\xi(p) = \chi_{[0,1]}(p)$. La distribución de la muestra es $X \sim Geo(p)$ luego $f_p(x) = p(1-p)^x$, luego la función de verosimilitud es \[ L(p|\textbf{x}) = f_n(\textbf{x}|p) = \prod_{i=1}^n f_p(x_i) = p^n (1-p)^{k} \] con $k = \sum_{i=1}^nx_i$. Como la función logaritmo es monótona creciente, maximizar esa función equivale a maximizar $\lambda_n(p|\textbf{x}) = \log L(p|\textbf{x}) = n \log(p) + k \log(1-p)$ y su derivada es $\frac{n}{p} +\frac{-k}{1-p}$ que es igual a cero si y solo si \[n(1-p)-kp=0 \iff p(k+n) = n \iff p = \frac{n}{k+n}\] luego el MLE es $p = \frac{n}{k+n}$.

Hallemos ahora con la fórmula de Bayes la distribución de $p|\textbf{x}$:\[ \xi(p|\textbf{x}) = \frac{\xi(p) f_n(\textbf{x}|p)}{f_n(\textbf{x})} = c(1-p)^kp^n \] con c una constante que no depende de $p$, luego $p|\textbf{x} \sim Beta(n-1,k-1)$.

Por tanto, el estimador bayesiano es $p_b = E(p|\textbf{x}) = \frac{n-1}{n+k-2}$.

Definamos ahora una función “estimador bayesiano cuadrático” análoga a la del apartado anterior:

  ebc<-function(k,n,L){     # función "estimador bayesiano cuadrático"
                          # si hay los mismos k la distribución es la misma
  pb<-(n+1)/(k+n+2)       # estimador bayesiano
  vp<-c()
  for(i in 1:L){
    u<-0
    while(u!=k){
      p<-runif(1,0,1)     # tomamos el parámetro a priori
      u<-rnbinom(1,n,p)   # sumamos muestras de geométricas supuesto el parámetro
    }                     # si la suma vale k me quedo con ese valor de p, si no me vale lo intento otra vez
    vp[i]<-p
  }                       # lo hago L veces
  print(n/(k+n))          # imprimo el MLE ( estimador de máxima verosimilitud)
  print(pb)               # imprimo el estimador bayesiano
  return(mean(vp))        # imprimo el promedio de los valores deseados
  }
  
  ebc(9,20,10000)

  ## [1] 0.6896552
  ## [1] 0.6774194

  ## [1] 0.6775582

Definamos ahora una función análoga al “error cuadrático medio gráfico” del apartado previo:

  ecmg<-function(k,n,L){     # error cuadrático medio gráfico
  pb<-(n+1)/(k+n+2)        # estimador bayesiano
  pm<-n/(k+n)              # MLE
  vp<-c()                  # idem construcción
  for(i in 1:L){
    u<-0
    while(u!=k){
      p<-runif(1,0,1)
      u<-rnbinom(1,n,p)
    }
    vp[i]<-p
  }
  puntos<-seq(0,1,0.01)               # una partición del intervalo (0,1) con norma 0.01
  ecm<-function(t)mean((vp-t)^2)      # función error cuadrático medio
  valores<-sapply(puntos,ecm)         # aplicamos la función a los puntos de la partición
  plot(puntos,valores,type='l',xlab='x',ylab='ecm') # esbozamos una curva con puntos en abscisas y valores en ordenadas
  segments(pm,0,pm,0.1,col='red')     # abscisa del MLE
  segments(pb,0,pb,0.1,col='green')   # abscisa el bayesiano
  print(ecm(pb))                      # el error cuadratico medio en el punto pb
  print((n+1)*(k+1)/((n+1+k+1)^2*(n+1+k+2)))        # la varianza de la distribución a posteriori
  print(var(vp))                      # la varianza del vector con los p deseados
  }
  
  ecmg(8,13,10000)

  ## [1] 0.009883572
  ## [1] 0.009924386
  ## [1] 0.009884257

En una distribución uniforme

Ejercicio 5

Sea ahora $X \sim Unif(0,\theta)$ con la distribución de $\theta$ de la forma $\theta \sim Unif(0,1)$. Luego la función de densidad a priori de $\theta$ es $\xi(\theta) = \chi_{[0,1]}(\theta)$.

Tenemos que $f_\theta (x) = \frac{1}{\theta} \chi_{[0,\theta]}(x)$, luego la función de verosimilitud es de la forma \[ L(\theta | \textbf{x}) = f_n(\textbf{x}|\theta) = \prod_{i=1}^n f_\theta (x_i) = \frac{1}{\theta^n}\chi_{[0,\theta]}(x_1) \cdot \chi_{[0,\theta]}(x_2) \cdots \chi_{[0,\theta]}(x_n) = \frac{1}{\theta^n}\chi_{[0,\theta]}(M) \] con $M = \max_{0\leq i \leq n}{\textbf{x}}$. Es claro que el MLE es el máximo $M$, ya que es el primer valor para el tenemos la certeza de que es el primer valor donde la función característica es no nula y la función característica es estrictamente monótona decreciente.

Empleando la fórmula de estimación de Bayes tenemos que la distribución a posteriori es de la forma\[ \xi(\theta|\textbf{x}) = \frac{\xi(\theta) f_n(\textbf{x}|\theta)}{f_n(\textbf{x})} = \frac{1}{\theta^nf_n(\textbf{x})} \chi_{[0,1]}(\theta) \chi_{[0,\theta]}(M) = \frac{1}{\theta^nf_n(\textbf{x})} \chi_{[M,1]}(\theta) \]Y como es una función de densidad conocemos que esta integral tiene que ir a 1, con esta idea hallemos la distribución marginal de la muestral: \[ 1 = \int_{0}^1 \frac{1}{\theta^nf_n(\textbf{x})} \chi_{[M,1]}(\theta) \, d\theta \iff f_n(\textbf{x}) = \int_{M}^1 \frac{1}{\theta^n} d\theta \, \]

analicemos lo que vale esta función para los distintos casos:

Caso $n=1$ \[ f_n(\textbf{x}) = \int_{M}^1 \frac{1}{\theta} d\theta \, =\log \left(\frac{1}{M}\right) \implies \xi(\theta|\textbf{x}) = \frac{1}{\theta^n f_n(\textbf{x})} \chi_{[M,1]}(\theta) = \frac{1}{\theta \log \left(\frac{1}{M}\right) } \chi_{[M,1]}(\theta) \]
Caso $n \geq 2$\[ f_n(\textbf{x}) = \int_{M}^1 \frac{1}{\theta^n} d\theta \, = \frac{\theta^{-n+1}}{-n+1} \Big|_{\theta=M}^{\theta=1} = \frac{1}{-n+1} -\frac{M^{-n+1}}{-n+1} = \frac{1}{n-1} \left(\frac{1}{M^{n-1}} - 1 \right) = \frac{1}{n-1}\frac{1-M^{n-1}}{M^{n-1}} \] \[ \implies \xi(\theta|\textbf{x}) = \frac{1}{\theta^n f_n(\textbf{x})} \chi_{[M,1]}(\theta) = \frac{1}{\theta^n \frac{1}{n-1}\frac{1-M^{n-1}}{M^{n-1}}} \chi_{[M,1]}(\theta) = \frac{(n-1)M^{n-1}}{ (1-M^{n-1}) \theta^n} \chi_{[M,1]}(\theta) \]

Hallemos ahora la esperanza de la distribución en distintos casos, es decir, el estimador bayesiano. Dividimos el cálculo de $E(\theta|\textbf{x})$en:

Caso $n=1$: \[E(\theta|\textbf{x}) = \frac{1}{\log\frac{1}{M}} \int_{M}^{1} \theta \frac{1}{\theta} \;d\theta = \frac{1-M}{\log\frac{1}{M}}\]
Caso $n=2$: \[ E(\theta|\textbf{x}) = \frac{M}{1-M} \int_{M}^{1} \theta \frac{1}{\theta^2} \;d\theta = \frac{M\log\left(\frac{1}{M}\right)}{1-M} \]
Caso $n>2$: \[ E(\theta|\textbf{x}) = \frac{(n-1)M^{n-1}}{1-M^{n-1}} \int_{M}^{1} \theta \frac{1}{\theta^n} \;d\theta = \frac{(n-1)M^{n-1}}{1-M^{n-1}} \left( \frac{1}{-n+2} \left( 1 - \frac{1}{M^{n-2}} \right)\right) = \frac{(n-1)}{(n-2)} M \frac{(1-M^{n-2})}{(1-M^{n-1})} \]

Comprobemos este resultado experiementalmente:

ebc<-function(M,n,L){            # función "estimador bayesiano cuadrático"
  vp<-c()
  for(i in 1:L){
    u<-0
    while(abs(max(u)-M)>0.01){
      p<-runif(1,0,1)           # tomamos el parámetro a priori
      u<-runif(n,0,p)           # sumamos n muestras uniformes supuesto el parámetro
    }                           # si el máximo de la muestra es aproximadamente el máximo deseado me lo quedo
    vp[i]<-p
  }
  if(n == 1) print( (1-M) / log(1/M))                          # estimador bayesiano cuando n=1
  else if(n == 2) print( M*log(1/M) / (1-M))                   # estimador bayesiano cuando n=2
  else print( (n-1)/(n-2) * M * (1-M^{n-2})/(1-M^{n-1}))       # estimador bayesiano cuando n>2
  print(M)                                                     # el MLE
  return(mean(vp))                                             # el promedio de los valores deseados
}
ebc(0.5,1,2000)

## [1] 0.7213475
## [1] 0.5

## [1] 0.7212587

ebc(0.7,2,2000)

## [1] 0.8322415
## [1] 0.7

## [1] 0.8303488

ebc(0.66,5,2000)

## [1] 0.7738371
## [1] 0.66

## [1] 0.7685618

Hallemos la esperanza de la función de penalización que será el error cuadrático medio, es decir, la varianza de la distribución a posteriori.

Caso $n=1$: \[ E(\theta^2|\textbf{x}) = \frac{1}{\log\frac{1}{M}} \int_{M}^{1} \theta^2 \frac{1}{\theta} \;d\theta = \frac{(1-M)^2}{2 \log\frac{1}{M}} \] \[\implies Var(\theta|\textbf{x}) = E(\theta^2|\textbf{x}) - \left(E(\theta|\textbf{x})\right)^2 = \frac{(1-M)^2}{2 \log\frac{1}{M}} - \left(\frac{1-M}{\log\frac{1}{M}}\right)^2 \]
Caso $n=2$: \[ E(\theta^2|\textbf{x}) = \frac{M}{1-M} \int_{M}^{1} \theta^2 \frac{1}{\theta^2} \;d\theta = M \] \[ \implies Var(\theta|\textbf{x}) = E(\theta^2|\textbf{x}) - \left(E(\theta|\textbf{x})\right)^2= M - \left(\frac{M\log\left(\frac{1}{M}\right)}{1-M}\right)^2 \]
Caso $n=3$: \[ E(\theta^2|\textbf{x}) = \frac{2M^{2}}{1-M^{2}} \int_{M}^{1} \theta^2 \frac{1}{\theta^3} \;d\theta = \frac{2M^{2}}{1-M^{2}} \log\left(\frac{1}{M}\right) \] \[ \implies Var(\theta|\textbf{x}) = E(\theta^2|\textbf{x}) - \left(E(\theta|\textbf{x})\right)^2 = \frac{2M^{2}}{1-M^{2}} \log\left(\frac{1}{M}\right) - \left( 2M \frac{(1-M)}{(1-M^{2})} \right)^2 \]
Caso $n>3$: \[ E(\theta^2|\textbf{x}) = \frac{(n-1)M^{n-1}}{1-M^{n-1}} \int_{M}^{1} \theta^2 \frac{1}{\theta^n} \;d\theta = \frac{n-1}{n-3}M^2 \frac{1-M^{n-3}}{1-M^{n-1}} \] \[ \implies Var(\theta|\textbf{x}) = E(\theta^2|\textbf{x}) - \left(E(\theta|\textbf{x})\right)^2 = \frac{n-1}{n-3}M^2 \frac{1-M^{n-3}}{1-M^{n-1}} - \left( \frac{(n-1)}{(n-2)} M \frac{(1-M^{n-2})}{(1-M^{n-1})} \right)^2 \]

Definamos ahora una función análoga al “error cuadrático medio gráfico” del apartado previo para comprobar los valores obtenidos:

ecmg<-function(M,n,L){     # error cuadrático medio gráfico
  if(n == 1) pb<-((1-M) / log(1/M))                        # estimador bayesiano cuando n=1
  else if(n == 2) pb<-(M*log(1/M) / (1-M))                 # estimador bayesiano cuando n=2
  else pb<- ((n-1)/(n-2) * M * (1-M^{n-2})/(1-M^{n-1}))    # estimador bayesiano cuando n>2
  pm<- M                                                   # MLE  
  vp<-c()                                                  # idem construcción
  for(i in 1:L){
    u<-0
    while(abs(max(u)-M)>0.01){
      p<-runif(1,0,1)           
      u<-runif(n,0,p)           
    }                           
    vp[i]<-p
  }
  puntos<-seq(0,1,0.01)               # una partición del intervalo (0,1) con norma 0.01
  ecm<-function(t)mean((vp-t)^2)      # función error cuadrático medio
  valores<-sapply(puntos,ecm)         # aplicamos la función a los puntos de la partición
  plot(puntos,valores,type='l',xlab='x',ylab='ecm') # esbozamos una curva con puntos en abscisas y valores en ordenadas
  segments(pm,0,pm,0.1,col='red')     # abscisa del MLE
  segments(pb,0,pb,0.1,col='green')   # abscisa el bayesiano
  print(ecm(pb))                      # el error cuadratico medio en el punto pb
  # la varianza de la distribución a posteriori
  if(n==1) print((1-M)^2/(2*log(1/M)) - ((1-M)/(log(1/M)))^2)        
  if(n==2) print(M - (M*log(1/M)/(1-M))^2)
  if(n==3) print((2*M^2/(1-M^2))*log(1/M) - (2*M*(1-M)/(1-M^2))^2)
  if(n>3) print((n-1)/(n-3)*M^2 * (1-M^(n-3))/(1-M^(n-1)) - ( (n-1)/(n-2)* M *(1-M^(n-2)/(1-M^(n-1))))^2)
  print(var(vp))                      # la varianza del vector con los p deseados
}
ecmg(0.7,2,1000)

## [1] 0.007500386
## [1] 0.007374026
## [1] 0.007487115

ecmg(0.8,3,1000)

## [1] 0.003156189
## [1] 0.003275837
## [1] 0.00315875

Ejercicio 6

El estimador bayesiano $p_{ba}$ que minimiza, en lugar del error cuadrático medio (es decir, con función de penalización $(\theta-a)^2$) el error absoluto medio (con función de penalización $|\theta-a|$) es la mediana.

Hallemos los distintos valores de la mediana:

Caso $n=1$: \[ \frac{1}{\log\frac{1}{M}} \int_{M}^{m} \frac{1}{\theta} \;d\theta = \frac{1}{2} \iff \frac{\log m +\log\frac{1}{M}}{\log\frac{1}{M}} = \frac{1}{2} \iff \frac{\log m}{\log\frac{1}{M}} + 1 = \frac{1}{2} \iff \log m= \frac{1}{2} \log M \iff m = \sqrt{M} \]
Caso $n \geq 2$: \[ \frac{(n-1)M^{n-1}}{ (1-M^{n-1})} \int_{M}^{m}\frac{1}{\theta^n} \; d\theta = \frac{1}{2} \iff \frac{(n-1)M^{n-1}}{ (1-M^{n-1})} \frac{1}{n-1} \left( \frac{1}{M^{n-1}} - \frac{1}{m^{n-1}} \right) = \frac{1}{2} \] \[ \iff \frac{1}{M^{n-1}} - \frac{1}{m^{n-1}} = \frac{1-M^{n-1}}{2M^{n-1}} \iff \frac{1}{m^{n-1}} = \frac{1+M^{n-1}}{2M^{n-1}} \iff m = M \left( \frac{2}{1+M^{n-1}} \right)^{\frac{1}{n-1}} \]

Veamos ahora el resultado con R

ebm<-function(M,n,L){            # función "estimador bayesiano que minimimiza el error absoluto medio"
  vp<-c()
  for(i in 1:L){
    u<-0
    while(abs(max(u)-M)>0.01){  # idem construcción
      p<-runif(1,0,1)           
      u<-runif(n,0,p)           
    }                           
    vp[i]<-p
  }
  if(n == 1) print( sqrt(M))                             # estimador bayesiano cuando n=1
  else if(n > 1) print( M *(2/(1+M^(n-1)))^(1/(n-1)) )   # estimador bayesiano cuando n>1                              
  return(median(vp))                                     # la mediana de los valores deseados
}                                          
ebm(0.5,1,2000)

## [1] 0.7071068

## [1] 0.708828

ebm(0.7,2,2000)

## [1] 0.8235294

## [1] 0.8250872

ebm(0.66,5,2000)

## [1] 0.7515152

## [1] 0.7544782

Práctica 3: Estimación Bayesiana

Esteban Chacobo, Víctor

2023-06-05