# P.d.f for Gamma posterior with Exponential data

I am trying to perform a simple exercise:

1. Sample $N$ points from $\text{Exponential}(\lambda=0.1)$
2. Assume a $\text{Gamma}(\alpha,\beta)$ prior for the parameter $\lambda$ above
3. Build a p.d.f for the posterior distribution, which is known to be $\text{Gamma}(\alpha+N,\beta+\sum_i^Nx_i)$ where $x_i$ is the sampled point.

Should I expect the resulting p.d.f to be peaked around the $\lambda = 0.1$, the parameter i used to generated the data? I understand that the initial parameters of the prior should have less and less impact on the posterior distribution with increasing number of sampled points.

The resulting histogram shows p.d.f that is far from being centered around $\lambda$. Is it a coding error or is my understanding flawed?

The python code goes:

import math
from numpy.random import gamma
from numpy.random import exponential
import matplotlib.pyplot as pp

N = 100
beta = 0.1
alpha = Shape = 0.1
Scale = Lambda = 1/beta

samples = exponential(scale=Scale, size = N)

dist = gamma(shape = alpha + N, scale = 1/beta + sum(samples), size = N)
count, bins, ignored = pp.hist(dist, 10, normed=True)

posterior_mean = (alpha+N)/(beta+sum(samples))

print("Posterior mean for the lambda parameter: {0}".format(posterior_mean))


In the above:

Shape $\equiv\alpha$

Scale $\equiv\frac{1}{\beta}\equiv\ \lambda$

The resulting histogram: Thanks!

• If $\lambda$ is known to be 0.1, how can you put a prior on it? – Greenparker Apr 2 '16 at 14:02
• I pretend I don't know the $\lambda$. I act as if the exponentially distributed dataset is given to me and i try to find the posterior distribution for the parameter. – ambushed Apr 2 '16 at 14:04
• I am not sure what coding language you have provided, but can you show us the output from the code, so we know exactly what is going on? – Greenparker Apr 2 '16 at 14:05
• In the code, it probably should be scale = 1/ (beta + sum(samples) ). You are missing parenthesis. – Greenparker Apr 2 '16 at 14:10
• Awesome, thanks! I wish i could mark your answer as the correct one as well. – ambushed Apr 2 '16 at 14:12

The Gamma distribution $\text{Gamma}(\alpha,\beta)$ has a mode at $\dfrac{\alpha-1}{\beta}$ and has a mean of $\dfrac{\alpha}{\beta}$
so if $N \gg \alpha$, $\sum_i^Nx_i \gg \beta$ and $\frac{1}{N} \sum_i^Nx_i \approx 10$
then $\text{Gamma}(\alpha+N,\beta+\sum_i^Nx_i)$ will have a mode and mean near $0.1$, so if you code is not producing this then you probably have a coding error. Greenparker has pointed out a possibility