# 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?

• I don't know how this affects the posterior but N should be discrete for the binomial and your prior for N is continuous? – einar Oct 9 '20 at 10:47
• Yes, N should be discrete, and I just used TruncatedNormal because I didn't know what else to use. This is literally all I know of my apple trees: N have been estimated previously to [6, 38, 53, 112, 140, 145, 264]. – Pål GD Oct 9 '20 at 11:52
• When I tried with a flat discrete prior, it over-estimated much more, especially if the upper bound was too high (e.g. 1000). – Pål GD Oct 9 '20 at 11:53
• I believe it over-estimates because there are two industry standard ways of estimating, and they both give lower numbers than my MCMC. – Pål GD Oct 9 '20 at 12:27
• Although you write you "end up with an estimate," I cannot find any description of how you make that estimate. What is your method?? Running a simulation cannot succeed because the logic is circular: in order to simulate you need to supply estimates in the first place. – whuber Oct 9 '20 at 13:59

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 $$$$q_{1:K} = (q_1, \ldots, q_K) .$$$$ The likelihood is $$$$p(q_{1:K}|N,\theta) = p(q_1|N,k)\prod_{k=2}^K p(q_k|q_{1:k-1},N,\theta) ,$$$$ 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: $$$$p(N,\theta|q_{1:K}) \propto p(q_{1:K}|N,\theta) .$$$$
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: $$$$\theta^{(r+1)} \sim p(\theta|q_{1:K},N^{(r)}) = \textsf{Beta}(\theta|a,b^{(r)})$$$$ 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 $$$$N^{(r+1)} = \begin{cases} N' & M \ge u \\ N^{(r)} & \text{otherwise} \end{cases} ,$$$$ where $$u \in \textsf{Uniform}(0,1)$$ and $$$$M = \frac{p(q_{1:K}|N',\theta^{(r+1)})}{p(q_{1:K}|N^{(r)},\theta^{(r+1)})} .$$$$ 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$$.
• Notational quibble: should it not be $p(q_k | N, q_1, \ldots, q_{k-1}, \theta) \sim Binomial(.)$? The full likelihood being p(q_{1:k}|N,\theta), with $q_i$s not independent but can be factorized $p(q_1 | N, \theta)p(q_2 | N, q_1, \theta),$ etc.? – einar Oct 14 '20 at 8:52