Poisson model

1 Exponential family

A family of distributions \(\mathcal F=\{p(y|\theta):\theta\in\Theta\subset\mathbb R^d\}\) is an \(\textit{exponential family}\) if its elements can be written as \[ p(y|\theta) = f(y)g(\theta) e^{\phi(\theta)^Tu(y)} \] where

  • \(f:\mathbb R\rightarrow \mathbb R\)
  • \(g:\mathbb R^d\rightarrow \mathbb R\)
  • \(\phi:\mathbb R^d\rightarrow \mathbb R^d\)
  • \(u:\mathbb R^d\rightarrow \mathbb R^d\)

are known functions. Prove that the Poisson distribution, with PMF

\[p(y|\theta) = \frac{\theta^{y}e^{-\theta}}{y!}, \quad \text{for } y=0,1,2,\ldots\]

with \(\theta>0\), is an exponential family distribution.

Sol

We can write the pdf as

\[p(y|\theta) = \frac{1}{y!} e^{y \log \theta -\theta}\] Thus, \(f(y)=\frac{1}{y!}\), \(g(\theta)=e^{-\theta}\), \(\phi(\theta)=\log(\theta)\), \(u(y)=y\).

2 Conjugate prior

Show that the Gamma prior \(\Gamma(\alpha,\beta)\) on \(\theta\) with pdf

\[\pi(\theta)= \frac{\beta^{\alpha}}{\Gamma(\alpha)}\theta^{\alpha-1}e^{-\beta\theta}, \quad \text{for } \theta > 0\]

with \(\alpha >0\) and \(\beta> 0\), is a conjugate prior for the Poisson model and find the parameters of the posterior distribution \(\pi(\theta|\textbf{y})\).

Sol

Conjugacy requires the posterior distribution \(\pi(\theta|\mathbf{y})\) and the prior distribution \(\pi(\theta)\) to be in the same family.

The posterior in our case is \[\theta|\mathbf{y} \sim Gamma\left(\sum_{i=1}^n y_i+\alpha, n+\beta\right)\]

3 Jeffreys’ prior

Compute the Jeffreys’ prior for the parameter \(\theta\) of a Poisson distribution. Recall that the Jeffreys’ prior in the univariate case is defined as \[\pi_J(\theta) \propto \sqrt{\mathcal I(\theta)}\] with \[\mathcal I(\theta) = - \mathbb{E}\left[\frac{d^2 \log p(\mathbf{y}|\theta)}{d \theta^2}\right]\] Is it a proper distribution? Compute the posterior distribution.

Sol

The Jeffreys’ prior is

\[\pi_J(\theta) \propto \theta^{-1/2}\]

which is an improper distribution. But the posterior distribution is proper and in particular

\[\pi(\theta|\mathbf{y}) \sim Gamma\left(\sum_{i=1}^n y_i+.5,n\right)\]

4 R Lab

Simulate three random samples of size 20, 100 and 1000 from a Poisson distribution with parameter \(\theta=.5\).

Assume a prior gamma distribution for the parameter \(\theta\) with the following four sets of hyper-parameters:

  • \(\alpha_1 = 1 \quad \beta_1 = .5\)
  • \(\alpha_2 = 1 \quad \beta_2 = 2\)
  • \(\alpha_3 = 1 \quad \beta_3 = 10\)
  • \(\alpha_4 = 2 \quad \beta_4 = 2\)
  1. plot the prior distributions Plot the four prior distributions on the same graph, then add the Jeffreys’ prior.
set.seed(123)
y20   <- rpois(20, .5)
y100  <- rpois(100, .5)
y1000 <- rpois(1000, .5)

a <- c(1, 1, 1, 2, .5)
b <- c(.5, 2, 10, 2, 0)

priordist <- function(theta, k){
    if(k==5){return(1/sqrt(theta))}
    return(dgamma(theta, a[k], b[k]))
}

xl <- c(0,3)
yl <- c(0,3)

plot(c(),c(),xlim=xl,ylim=yl,main="Prior",xlab="",ylab="")
for(k in 1:5) curve(priordist(x,k), add=T, col=k)
legend(2.5, 2.5,
       legend=c("Gamma1", "Gamma2",
                "Gamma3", "Gamma4",
                "Jeffreys'"),
       col=c(1:5), lty=1, cex=0.8)

  1. Compare the posterior distributions For each simulated data set, plot on the same graph the five posterior distributions.
postdist <- function(theta,y,k){
  return(dgamma(theta,a[k]+sum(y),b[k]+length(y)))
}

plot(c(),c(),xlim=xl,ylim=yl,main="n=20",xlab="",ylab="")
for(k in 1:5) curve(postdist(x,y20,k), add=T, col=k)

plot(c(),c(),xlim=xl,ylim=yl,main="n=100",xlab="",ylab="")
for(k in 1:5) curve(postdist(x,y100,k), add=T, col=k)

plot(c(),c(),xlim=xl,ylim=yl,main="n=1000",xlab="",ylab="")
for(k in 1:5) curve(postdist(x,y1000,k), add=T, col=k)

  1. Find the posterior expectations and the maxima a posteriori (MAP) of \(\theta\).
## posterior expectations
(a+sum(y20))/(b+20)
## [1] 0.5853659 0.5454545 0.4000000 0.5909091 0.5750000
(a+sum(y100))/(b+100)
## [1] 0.4975124 0.4901961 0.4545455 0.5000000 0.4950000
(a+sum(y1000))/(b+1000)
## [1] 0.4797601 0.4790419 0.4752475 0.4800399 0.4795000
##MAP
(a+sum(y20)-1)/(b+20)
## [1] 0.5365854 0.5000000 0.3666667 0.5454545 0.5250000
(a+sum(y100)-1)/(b+100)
## [1] 0.4875622 0.4803922 0.4454545 0.4901961 0.4850000
(a+sum(y1000)-1)/(b+1000)
## [1] 0.4787606 0.4780439 0.4742574 0.4790419 0.4785000
  1. Construct the 95% credible intervals, equi-tails and highest posterior density (HPD), for the three data sets using the Jeffreys’ prior.
hpd<-function(y,p){
dy<-density(y)
md<-dy$x[dy$y==max(dy$y)]
py<-dy$y/sum(dy$y)
pys<--sort(-py)
ct<-min(pys[cumsum(pys)< p])
list(hpdr=range(dy$x[py>=ct]),mode=md)
}

hpd(rgamma(1000,-.5+sum(y20),20),.95)
## $hpdr
## [1] 0.2171316 0.8405557
## 
## $mode
## [1] 0.4906497
hpd(rgamma(1000,-.5+sum(y100),100),.95)
## $hpdr
## [1] 0.3434509 0.6333469
## 
## $mode
## [1] 0.4596539
hpd(rgamma(1000,-.5+sum(y1000),1000),.95)
## $hpdr
## [1] 0.4365653 0.5209245
## 
## $mode
## [1] 0.4770072
qgamma(c(.025,.975),-.5+sum(y20),20)
## [1] 0.2570724 0.8869719
qgamma(c(.025,.975),-.5+sum(y100),100)
## [1] 0.3582076 0.6307072
qgamma(c(.025,.975),-.5+sum(y1000),1000)
## [1] 0.4365814 0.5223127