# Gamma-Poisson conjugate prior, posterior exploding?

I've been looking for simple code that can model ad clicks per day. Notionally, gamma-poisson would be a good conjugate prior. However, I'm finding that for slightly large daily click rate values, the denominator, (n-1)! explodes.

I'd like to know:

(A) Am I doing it right?

(B) Is there a better conjugate prior to choose?

My code:


# section 1: imports
from scipy.stats import poisson
from scipy.stats import gamma as gamma_dist
import matplotlib.pyplot as plt
from math import factorial
from math import gamma
import scipy

# section 2: class definition
class Gamma():

def __init__(self,prior_mu,cutoff):
self.a = prior_mu
self.b = 1
self.span = np.linspace(0,cutoff,100)

def eval_density(self,a,x,b):
num = b**(a)*x**(a-1)*np.exp(-x*(b+1))
denom = gamma(a)
denom = scipy.special.gamma(a)
return num/denom

def update(self,batch):
self.a += batch
self.b += 1

def plot(self):
density = [self.eval_density(self.a,x,self.b) for x in self.span]
plt.plot(self.span,density)

# section 3: data generation and experiment definition
data = gamma_dist.rvs(a=4,scale = 1,size = 50)
def experiment_1(var,data,prop=0.2):
window = (len(data)*prop)
for idx,obs in enumerate(data):
try:
if idx%window==0:
var.plot()
var.update(obs)
except:
return f"overflow at {idx}"

# section 4: trial run
g = Gamma(cutoff=15,prior_mu=2)
experiment_1(var=g,data=data)


And the resulting plot:

As you can see from code, my prior belief was that the rate was 2 clicks per day. (In truth this is simulated data and the actual rate is 4.) The plot does slowly converge, however, the peak shrinks quite a bit and isn't necessary tightening the variance.

I've used similar code for a Beta-Binomial conjugate prior and the results were night and day different. In the beta case, the peaks increased and became tighter with more data. In the gamma case, the peaks reduced and ultimately the code crashed after 40 of 50 iterations because the denominator exploded.

Feels like I'm doing it wrong.

A gamma prior with shape parameter $$\alpha_0 = 2000$$ and rate parameter $$\kappa_0 = 1),$$ for the Poisson mean $$\lambda$$ has $$E(\lambda) = 2000$$ and $$P(\lambda < 2100) \approx 0.99.$$ Maybe that is a reasonable prior for a click rate "about" $$2000,$$ but not likely more than 2100.

pgamma(2100, 2000, 1)
[1] 0.9863525


Subsequently, if data over $$n=20$$ days shows $$t=42\,000$$ clicks, then considering $$t = \sum_{i=1}^{20} x_i,$$ where $$x_i \stackrel{iid}{\sim} \mathsf{Pois}(\lambda),$$ the likelihood function is $$\lambda^t\,e^{-n\lambda}.$$

Thus, the posterior density is of the form $$p(\lambda|x) \propto \lambda^{\alpha_0-1}e^{\kappa_0\lambda} \times \lambda^t\,e^{-n\lambda} = \lambda^{\alpha_0+t-1}\,e^{-(\kappa_0+n)\lambda},$$ which is the kernel (density without norming constant) of $$\mathsf{Gamma}(\alpha_n, \kappa_n),$$ where $$\alpha_n = \alpha_0+t,\,\kappa_n=\kappa_0+n).$$

Thus, for our example, the posterior distribution is $$\mathsf{Gamma}(44000, 21),$$ the posterior mean is $$\alpha_n/\kappa_n \approx 2095$$ and a 95% Bayesian probability interval is $$(2075,2115).$$

a.n = 44000;  k.n = 21
qgamma(c(.025,.975), a.n, k.n)
[1] 2075.706 2114.861


If we were doing frequentist inference, then a Wald 95% confidence interval (reasonably useful because of the large numbers of counts) would be about $$(2080, 2120).$$

pm = c(-1,1); (42000 + pm*1.96*sqrt(42000))/20
[1] 2079.916 2120.084


If you want to choose a prior distribution that has less influence on the posterior, then choose both $$\alpha_0$$ and $$\beta_0$$ to be very small. A 95% Bayesian probability interval with a noninformative prior (say, using $$\alpha_0 = \kappa_0 = 0.01),$$ would be $$(2079, 2119).$$

a.n = 42000.01;  k.n = 20.01
qgamma(c(.025,.975), a.n, k.n)
[1] 2078.925 2119.072

• This might be a novice question but here goes- if we restrict our attention to just the numerator (not scaled by the normalizing constant) and it is proportional (but not equal to) the posterior -- do we have enough information to compute the expectation and variance? I think from your solution, the answer is yes.. But I feel a bit out of my depth. Jul 9, 2020 at 13:31
• Once can find posterior mean $(\approx 2100)$ and variance $(\approx 105).$ Sep 6, 2020 at 19:21