# posterior distribution of a Poisson mixture model

This is a Poisson-gamma model with mixture prior, thus mixture posterior. I am having some trouble finding the posterior weightings.

I have the prior weightings $$p_1=1/3$$; $$p_2=2/3$$. The 2 component priors are $$\theta_1$$~$$\Gamma_{220,10}$$, and $$\theta_2$$~$$\Gamma_{110,14}$$. I am given 7 data points with Poisson counts for each $$\{y_1,...,y_7\}$$. So $$y_1,...,y_7|\theta_j$$~ $$Poisson_{\theta_j}$$;j can be 1 or 2. $$\Sigma y_i= 189$$.

I have tried to derive the analytic form of the posterior, which looks like this:

The code I have done so far. The term 'gamma(b1+n-y)' is causing the troubles:

#a. mixture prior:
N=1000
a1<-200; b1<-10
a2= 110; b2=4
a<-c(a1,a2); b<-c(b1,b2)
# the prior weightings are (1/3, 2/3):
gamma=1/3

# MC simulation of mixture prior:
lambda<-seq(0,1,0.01)
g<-rbinom(N,1,gamma)+1  # generate random bernoulli rvs+1 = {1,2} according to prior weightings
gprior<-rgamma(N,a[g],b[g]) # g being 1 or 2 indicates the parameter to use in rbeta

# plot prior:
par(mfrow=c(1,1))
plot(density(gprior),xlab="lambda",lwd=2,lty=2,col="gray",main="Mixture of conjugate priors",cex.main=0.8,cex.axis=0.8)

#b. posterior:
data= c(28,26,26,19,22,38,30)
n= length(y)
y= sum(data)

# The posterior normalising constant (inverse):
k= (10^220)*gamma(409)/(3*gamma(220)*17^409)+ 2*4^110*gamma(299)/(3*gamma(110)*11^299)
gamma.star= (1/k)*10^220*gamma(409)/(3*gamma(220)*17^409)

# use the log-gamma to rescale:
k= log(10^220)*lgamma(409)/(3*lgamma(220)*log(17^409))+ 2*log(4^110)*lgamma(299)/(3*lgamma(110)*log(11^299))
gamma.star= (1/k)*log(10^220)*lgamma(409)/(3*lgamma(220)*log(17^409))


• The posterior weights are not correct. – Xi'an Aug 28 '20 at 4:35
• @Xi'an Yes thank you for pointing out. I realise that those weights are for the binomial sampling model. However as I tried straight derivation, the huge gamma functions and exponents are giving me troubles. Rescaling by log() and lgamma() does not entirely solve the problem as log(17^409) is still Inf-valued. What can I do about that? – siegfried Aug 28 '20 at 8:13

When computing the posterior weights the normalising constant is not needed: \begin{align*} \gamma^* &\propto \frac{1}{3}\frac{10^{220}\Gamma(409)}{\Gamma(220)17^{409}}\\ 1-\gamma^* &\propto \frac{2}{3}\frac{4^{110}\Gamma(299)}{\Gamma(110)11^{299}} \end{align*} means that $$\gamma^*=\dfrac{\frac{1}{3}\frac{10^{220}\Gamma(409)}{\Gamma(220)17^{409}}}{\frac{1}{3}\frac{10^{220}\Gamma(409)}{\Gamma(220)17^{409}}+\frac{2}{3}\frac{4^{110}\Gamma(299)}{\Gamma(110)11^{299}}}$$ can be simplified as $$\gamma^*=\dfrac{1}{1+2\frac{4^{110}\Gamma(299)}{\Gamma(110)11^{299}}\frac{\Gamma(220)17^{409}}{10^{220}\Gamma(409)}}$$ and the term in the denominator remains manageable:
> lgamma(299)+lgamma(220)-lgamma(110)-lgamma(409)+

• Yes,$$\gamma^*=\dfrac{1}{1+2\exp\{1.915\}}$$with the data provided in the question – Xi'an Sep 13 '20 at 16:11