0
$\begingroup$

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: enter image description here

Please help me with the posterior weightings. Thanks in advance.

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))



$\endgroup$
  • $\begingroup$ The posterior weights are not correct. $\endgroup$ – Xi'an Aug 28 '20 at 4:35
  • $\begingroup$ @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? $\endgroup$ – siegfried Aug 28 '20 at 8:13
1
$\begingroup$

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)+
  110*log(4)+409*log(17)-299*log(11)-220*log(10)
[1] 1.915338
$\endgroup$
  • $\begingroup$ Hi I would just like to confirm that I will need to take exp to 1.915338 before calculating the gamma* term? My lecturer changed the data magnitude in the previous problem so I did not use your method. But here I am with some huge real data following a bimodal Poisson model... $\endgroup$ – siegfried Sep 13 '20 at 12:27
  • 1
    $\begingroup$ Yes,$$\gamma^*=\dfrac{1}{1+2\exp\{1.915\}}$$with the data provided in the question $\endgroup$ – Xi'an Sep 13 '20 at 16:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.