Estimating the number of apples in an apple tree using MCMC I'm trying to estimate the number of apples in an apple tree by
repeatedly kicking the tree and counting how many apples fall down.
This process, I believe, is called removal sampling.
The only assumption I'm making is that there is a constant probability
$p$ that an apple falls down when I kick the tree.
Given that, e.g., [100, 10, 1, 0] apples fall down, I'm prone to
believe that $p \sim 0.9$ and that $N$, the total number of apples, is
111.
I'm trying to estimate $N$ using MCMC, but I seem to struggle to get it
right.
First, I've picked $N$ to be a uniform discrete prior in $[0, 1000]$ and $p$ is chosen to be a flat uniform over $[0, 1]$.
But, suppose that I remove, e.g., [19, 17, 13, 1, 1], I used
$$
\begin{align}
N   &= \text{DiscreteUniform}(0, 1000)\\
p   &= \text{Uniform}(0, 1)\\
q_1 &= \text{Binomial}(N, p, \text{observed}=19)\\
q_2 &= \text{Binomial}(N - q_1, p, \text{observed}=17)\\
q_3 &= \text{Binomial}(N - (q_1 + q_2), p, \text{observed}=13)\\
q_4 &= \text{Binomial}(N - (\sum_j^{3} q_j), p, \text{observed}=1)\\
q_5 &= \text{Binomial}(N - (\sum_j^{4} q_j), p, \text{observed}=1)
\end{align}
$$
With this model, I end up with an estimate on $N \sim 54.5$.  A different model I tried consistently get $N \leq 53$.
My question is the following:
Are my $q_i$ correctly modeled, or are they possibly the reason my model is over-estimating?
Do you see any obvious mistakes?


 A: Your results look reasonable given your model and your other assumptions. I can't speak to whether the model (and the assumptions) are themselves reasonable.
I'm going to change the notation a bit because I like to use the "$p$" to denote a probability density (or mass) function. So I'll use $\theta$ as the probability of "success" instead.
The observations are given by
\begin{equation}
q_{1:K} = (q_1, \ldots, q_K) .
\end{equation}
The likelihood is
\begin{equation}
p(q_{1:K}|N,\theta) = p(q_1|N,k)\prod_{k=2}^K p(q_k|q_{1:k-1},N,\theta) ,
\end{equation}
where
\begin{align}
p(q_1|N,\theta) &= \textsf{Binomial}(q_1|N,\theta) \\
p(q_k|q_{1:k-1},N,\theta) &= \textsf{Binomial}\left(q_k\Big|N - \sum_{j=1}^{k-1} q_j,\theta\right) .
\end{align}
The prior for the latent variables $N$ and $\theta$ is flat. Therefore, the posterior is proportional to the likelihood:
\begin{equation}
p(N,\theta|q_{1:K}) \propto p(q_{1:K}|N,\theta) . 
\end{equation}
We can apply the Metropolis-within-Gibbs sampler. Let $(\theta^{(r)},N^{(r)})$ denote the current state of the chain. The full-conditional posterior distribution for $\theta$ delivers the following:
\begin{equation}
\theta^{(r+1)} \sim p(\theta|q_{1:K},N^{(r)}) = \textsf{Beta}(\theta|a,b^{(r)}) 
\end{equation}
where
\begin{align}
a &= 1 + \sum_{k=1}^K q_k \\
b^{(r)} &= 1 + K\,N^{(r)} - \sum_{k=1}^K (K-k+1)\,q_k .
\end{align}
To draw $N^{(r+1)}$ we can take a Metropolis step, using a symmetric uniform discrete proposal. Let $N' = N^{(r)} + \delta$ where $\delta \sim \textsf{Uniform}(-3,3)$ for example. Then
\begin{equation}
N^{(r+1)} = \begin{cases}
N' & M \ge u \\
N^{(r)} & \text{otherwise}
\end{cases} ,
\end{equation}
where $u \in \textsf{Uniform}(0,1)$ and
\begin{equation}
M = \frac{p(q_{1:K}|N',\theta^{(r+1)})}{p(q_{1:K}|N^{(r)},\theta^{(r+1)})} .
\end{equation}
Given the draws $\{(\theta^{(r)},N^{(r)})\}_{r=1}^R$ one can produce figures that are similar to those in the question (except that the distribution for $N$ should be discrete). The posterior mode for $N$ is 53 and the posterior mean is about 55. There's about a 10% chance that $N\ge 60$.
