# Elementary MCMC pseudocode

My aim relates to the project described here, but I've tried to make this question self-contained.

I'm trying to write the MCMC pseudocode for the following inference problem:

Given two observed independent outcomes $X$ and $Y$, what are the posterior densities of probabilities $p_{x}$ and $p_{y}$ of obtaining $X$ and $Y$, assuming $X$ ~ $\text {Bin}(N, p_{x})$ and $Y$ ~ $\text{Bin}(N, p_{y})$? We do not know or observe $N$. We have priors for $N$, $p_{x}$, and $p_{y}$. We will assume that $N$ ~ $\text{Poisson}(\lambda)$, and $p_{x}$ and $p_{y}$ will have either $\text{Beta}$ or $\text{Uniform}$ priors, chosen by the user. For the $\text{Beta}$ priors, the mean and variance will be specified. For the $\text{Uniform}$ priors, the upper and lower limits will be specified.

My background is obviously not statistics, and I'm having a hard time patching gaps in my knowledge from the voluminous MCMC literature. Rather than burden commenters with explaining a mountain of Bayesian inference to me, I'd appreciate it if people could point out the errors in my pseudocode and understanding, including instances in which I might be misusing terms. I'll focus my study on those areas. I welcome suggested reading that is specific (I've seen but not completely read some of the suggestions here).

I will be coding this by hand from scratch, except for simple math functions.

Ultimately, I'll make sure there's complete and accurate pseudocode associated with this post.

Proposed pseudocode

1. Declare an initial value of $\theta_{0} = (N,p_{x},p_{y})$ and set equal to $\theta_{\text{current}}$.
2. Calculate the posterior probability of $\theta_{\text{current}}$:

• First, calculate the log-likelihood $\textit{L}_\text{current}$ of $\theta_{\text{current}}$. Here, the log-likelihood $\textit{L}_\text{current}$ = $L_{x} + L_{y}$, where $L_{x}$ $=$ $\text{log}{N \choose X} + X$ $\text{log } p_{x} + (N-X) \text{ log }$ $(1-p_{x})$ and $L_{y}$ is analogous.

• Next, obtain the posterior by multiplying the log-likelihood $\textit{L}_\text{current}$ by the prior, $\text{Prior}(\theta_\text{current})$. The prior of $\theta_\text{current}$ is the product of the priors of each of the parameters in $\theta$. In this case, if the user has chosen $\text{Beta}$ priors for $p_{x}$ and $p_{y}$, the respective priors are simply the $\text{Beta}$ pdf evaluated at $p_{x}$ or $p_{y}$. The prior for $N$ is the pmf of the Poisson function at $N$.

3. Draw $\theta_{\text{candidate}}$ from the proposal distribution. What exactly is this distribution for each parameter in $\theta$? I've seen normal distributions used, but here $N$ and $p_{x}$ clearly must be drawn from different distributions.
4. Calculate the posterior of $\theta_{\text{candidate}}$ following step 2.
5. Evaluate whether to accept $\theta_{\text{candidate}}$ by comparing the posteriors, and update $\theta_{\text{current}}$ if justified. I will use the Metropolis algorithm (not shown for brevity).
6. Repeat steps 3-5 until convergence or reasonable stopping point.

(tl;dr) Major questions

• In step 2, am I missing shortcuts related to the conjugate priors? I still don't really understand what they're good for or if I'm calculating the prior correctly.

• How are jump distributions determined?

-
The first thing I notice, though this might be a typo is: "The prior for $N$ is the pdf of the exponential function at $N$." Do you mean poisson? Also, as you assumed an exponential prior on $\lambda$, you need to add that as an additional variable. – Cam.Davidson.Pilon Jan 10 '13 at 20:27
Secondly, have you looked at Gibbs sampling? This is likely much easier for steps 3-5. – Cam.Davidson.Pilon Jan 10 '13 at 20:40
@Cam.Davidson.Pilon: Quite right about the prior for $N$. Typo. Fixing. It's $\lambda$, not $N$, that's being estimated. I knew more about Gibbs sampling 4-5 y ago and can look into it again... but if there's a decent way to work from what's here, I'd like to try it first. – neophyte Jan 11 '13 at 0:27
I commented too soon. It's not clear to me why both $N$ and $\lambda$ would be estimated (they're both the expected number of trials). – neophyte Jan 11 '13 at 0:37

• For the proposal distribution, just try a random walk. If a proposal is outside the support of $\theta$, then reject it. (Think of it as having a prior value of zero.)