Bayesian Aproach: Infering the N and $\theta$ values from a binomial distribution

Question

I am doing a homework about infering the N value of a binomial distribution for my Bayesian Statistics Course and I have seen a paper in Biometrika magazine published in 1988 for doing so. The question is how to simulate the N and $\theta$ values based on their respective marginal posteriors. The properties about the model are explained below:

$$X\sim Binom(N,\theta)$$ The priori distribution of N:$$ N\sim Poisson(\mu)$$ $$\lambda = \mu\theta$$(The logic here is that it may be easier to formulate a prior considering the unconditional expectation of the observations, rather than the mean of the unobserved N).
The $\lambda$ and $\theta$ a prior are distributed as:$$\lambda \sim Gamma(\kappa_{1},\kappa_{2})$$ $$\theta\sim Beta(\alpha,\beta)$$

After doing the procedings found on the paper which this link is below, the N marginal posterior is the following one: $p(N|x_{1},x_{2},..,x_{n})\propto(N!)^{-1}\Gamma(N+\kappa_1)\{\prod_{i=1}^5\dbinom{N}{x_i}\} \int_0^1 \! \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta^{-N+S+\alpha-1}(1-\theta)^{nN-S+\beta-1}(\theta^{-1}+\kappa_{2})^{(N+\kappa_{1})} \, \mathrm{d}\theta $.

How can create a random variable generator from that marginal probability function?

Then, I suppose that this integral:$ \int_0^1 \! \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta^{-N+S+\alpha-1}(1-\theta)^{nN-S+\beta-1}(\theta^{-1}+\kappa_{2})^{(N+\kappa_{1})} \, \mathrm{d}\theta $ represents part of the marginal posterior of $\theta$ but I don´t know if this is true, My second question is that: Which could be the marginal distribution function for $\theta$? and How can I create a random variable function for getting $\theta$ values?.

I attached the paper I 've working with, the only thing I added is a Beta prior distributión for $\theta$.

http://ftp.stat.washington.edu/raftery/Research/PDF/bka1988.pdf

Answer 1

The hierarchical model seems to be $$ \begin{eqnarray} X_i\mid N=n,\Theta=\theta,\Lambda=\lambda &\sim& \mathrm{Bin}(n,\theta) \qquad\qquad\qquad i=1,\dots,m \\ N\mid\Theta=\theta,\Lambda=\lambda &\sim& \mathrm{Poisson}(\lambda/\theta) \\ \Theta &\sim& \mathrm{Beta}(\alpha,\beta) \\ \Lambda &\sim& \mathrm{Gamma}(\kappa_1,\kappa_2) \, , \\ \end{eqnarray} $$ in which $\alpha,\beta,\kappa_1$, and $\kappa_2$ are known. Defining $X=(X_1,\dots,X_m)$ and $x=(x_1,\dots,x_m)$, and using Bayes' Theorem, we have $$ \begin{eqnarray} f_{N,\Theta,\Lambda\mid X}(n,\theta,\lambda\mid x) &\propto& f_{X\mid N,\Theta,\Lambda}(x\mid n,\theta,\lambda)\; f_{N,\Theta,\Lambda}(n,\theta,\lambda) \\ &=& \left(\prod_{i=1}^m f_{X_i\mid N,\Theta,\Lambda}(x_i\mid n,\theta,\lambda)\right) \; f_{N\mid\Theta,\Lambda}(n\mid\theta,\lambda) \; f_\Theta(\theta) \; f_{\Lambda}(\lambda) \, .\\ \end{eqnarray} $$ Therefore, $$ \begin{eqnarray} f_{N\mid X}(n\mid x) &=& \int_0^1 \int_0^\infty f_{N,\Theta,\Lambda\mid X}(n,\theta,\lambda\mid x)\,d\lambda\,d\theta \\ &\propto& \frac{1}{n!} \left( \prod_{i=1}^m {n\choose x_i}\right) \int_0^1 \int_0^\infty \theta^{t-n+\alpha -1} (1-\theta)^{mn-t+\beta-1} \\ && \qquad\qquad\qquad\qquad\qquad\times\lambda^{n+\kappa_1-1} \exp\left(-\left(\frac{1+\kappa_2\theta}{\theta}\right)\lambda\right)\,d\lambda\,d\theta \, , \end{eqnarray} $$ in which $t=\sum_{i=1}^m x_i$. The integration in $\lambda$ is easy. Just identify the kernel of a $$ \mathrm{Gamma}\left( n+\kappa_1, \frac{1+\kappa_2\theta}{\theta} \right) $$ distribution, and adjust the normalization constant. Doing this, we have $$ f_{N\mid X}(n\mid x) \propto \frac{\Gamma(n+\kappa_1)}{n!} \left( \prod_{i=1}^m {n\choose x_i}\right) \int_0^1 \frac{\theta^{t+\alpha-mn-1}(1-\theta)^{mn-t+\beta-1}}{(1+\kappa_2\theta)^{n+\kappa_1}} \,d\theta \, . $$ Of course, the last expression holds for $n\geq \max\{x_1,\dots,x_m\}$. If this is not clear, write the indicators $I_{\{0,\dots,n\}}(x_i)$ in the binomials explicitly.

If we use an improper prior $f_\Lambda(\lambda)\propto 1/\lambda$, the integrations are easy. Using an uniform prior $(\alpha=\beta=1)$ for $\Theta$, we find the closed form $$ f_{N\mid X}(n\mid x) \propto \frac{n\,\Gamma(t+1)\Gamma(mn-t+1)}{\Gamma(mn+2)} . $$

How do we normalize this? What is the Bayes estimate with quadratic loss $\mathrm{E}[N\mid X=x]$? Can we find analytically at least the posterior mode? Can we sample from this distribution to find an estimate and a credible set?

I didn't check the algebra. Please, do it.