Ejercicio Bayes

Ejemplo con distribución Poisson

Suponga que tenemos una muestra aleatoria $X_1,...,X_{10}$ tal que $X_j \sim Poisson(\theta)$. Recordemos que

1. Datos

set.seed(10000)
n<-10
X<-rpois(n, 4) #suponemos que theta = 4.
X

##  [1] 4 4 3 4 3 3 5 5 4 1

mean(X)

## [1] 3.6

Tenemos que la verosimilitud es:

$\mathcal{L}(\theta|x) = \frac{ \theta^{\sum_{j=1}^{n}x_j} e^{-n\theta} }{\prod_{j=1}^{n} x_{j}! }$

2. Distribución previa no informativa.

$$\theta \sim Unif(0,B)$$ donde $B$ va a representar un valor arbitrario muy grande.

La distribución a posteriori es:

$$\pi(\theta| x) \propto \mathcal{L}(\theta |x)\pi(\theta)$$

$$\Rightarrow \pi(\theta| x) \propto \frac{ \theta^{\sum x_j} e^{-n\theta} }{\prod x_{j}! } \cdot \frac{1}{B} \propto \theta^{n\bar{x}} e^{-n\theta}$$ Esto indica que la distribución posterior tiene nucleo de una distribución gamma:

$$\theta|x \sim Gamma(n\bar{x} + 1, n)$$

B<-100

alfa1<- n*mean(X)+1
beta1<- n

theta <- seq(0,10,by=0.1)
pi.theta <- dgamma(x=theta,shape=alfa1,rate=beta1)
plot(theta,pi.theta,type="l")
abline(v=mean(X),col=2)
abline(h=1/B,col=3)

La media posterior usando una previa no informativa.

alfa1/beta1

## [1] 3.7

3. Distribución previa informativa.

Suponga que conocemos de antemano información respecto a $\theta$,

$$\pi(\theta) = \frac{ \beta^{\alpha} \theta^{\alpha - 1} e^{-\beta \theta} }{ \Gamma(\alpha) }$$ $$\pi(\theta|x) = \frac{\mathcal{L}(\theta |x)\pi(\theta)}{\int_{0}^{+\infty} \mathcal{L}(\theta|x)\pi(\theta)} \propto \mathcal{L}(\theta |x)\pi(\theta) = \frac{ \theta^{n\bar{x}} e^{-n\theta} }{\prod x_{j}! } \cdot \frac{ \beta^{\alpha} \theta^{\alpha - 1} e^{-\beta \theta} }{ \Gamma(\alpha) }$$ $$\propto \theta^{n\bar{x}} e^{-n\theta} \cdot \theta^{\alpha - 1} e^{-\beta \theta} = \theta^{n\bar{x} + \alpha - 1} e^{-\theta(n + \beta)}$$

Por lo tanto concluimos que $\theta|x \sim Gamma(n\bar{x} + \alpha , n + \beta)$.

alfa0 <- 4
beta0 <- 2
  
alfa1<- n*mean(X)+alfa0
beta1<- n+beta0

theta <- seq(0,10,by=0.1)
pi.theta.prior <- dgamma(x=theta,shape=alfa0,rate=beta0)
pi.theta <- dgamma(x=theta,shape=alfa1,rate=beta1)
plot(theta,pi.theta.prior,type="l", ylim=c(0,1))
points(theta,pi.theta,type="l",col=2)
abline(v=mean(X),col=3)

La media posterior usando una previa informativa (depende de la información a priori).

alfa1/beta1

## [1] 3.333333

Recuerde que:

$$E(\theta|x) = \frac{n\bar{x} + \beta\left( \frac{\alpha}{\beta}\right)}{n + \beta}$$