Maximum likelihood estimator of $p$ for the binomial (truncated) distribution

Question

This question is from the book "Introduction To The Theory Of Statistics, Mood Alexander", Chapter 7, Problem 15. In genetic investigations one frequently samples from a binomial distribution, except that observations of $x=0$ are impossible; so, in fact, the sampling is from the conditional (truncated) distribution

\begin{equation*} \dbinom{m}{x} \dfrac {p^x(1-p)^{m-x}}{1-(1-p)^m}\mathbf{I}_{\{1,2,\cdots,m\}}(x) \end{equation*}

Find the maximum-likelihood estimator of $p$ in the case $m=2$ for samples of size $n$.

So, i did the calculations by find the zero of the first derivative from the likelihood function, then i got:

\begin{align*} L(p|x_{1}, x_{2}, \cdots, x_{n}) & = \displaystyle\prod_{i=1}^{n} \dbinom{2} {x_{i}}\dfrac{p^{x_{i}}(1-p)^{2-x_{i}}}{1-(1-p)^{2}} \mathbf{I}_{\{1,2\}}(x_{i}) \\ \mbox{log}\hspace{1ex}L(p|x_{1}, x_{2}, \cdots, x_{n}) & = \mbox{log} \left[\displaystyle\prod_{i=1}^{n} \dbinom{2}{x_{i}}\right] - n\mbox{log}[p(2-p)] + \displaystyle\sum_{i=1}^{n}x_{i}\mbox{log} (p)\\ &\qquad\qquad + \left( 2n - \displaystyle\sum_{i=1}^{n}x_{i} \right)\mbox{log}(1-p)\\ \dfrac{d}{dp}\mbox{log}\hspace{1ex}L(p|x_{1}, x_{2}, \cdots, x_{n}) & = - n \cdot\dfrac{(2-2p)}{[p(2-p)]} + \displaystyle\sum_{i=1}^{n}x_{i}\cdot\dfrac{1}{p}\\ &\qquad\qquad - \left( 2n - \displaystyle\sum_{i=1}^{n}x_{i} \right)\cdot\dfrac{1}{1-p}\\ & = \dfrac{- 2n + 2\displaystyle\sum_{i=1}^{n}x_{i} - p \displaystyle\sum_{i=1}^{n}x_{i}}{p(1-p)(2-p)} \\ \end{align*}

Now,

\begin{align*} \dfrac{d}{dp}\mbox{log}\hspace{1ex}L(p|x_{1}, x_{2}, \cdots, x_{n}) & = 0 \\ \dfrac{- 2n + 2\displaystyle\sum_{i=1}^{n}x_{i} - p \displaystyle\sum_{i=1}^{n}x_{i}}{p(1-p)(2-p)} & = 0 \\ \hat{p} & = 2\left( 1 - \dfrac{n}{\displaystyle\sum_{i=1}^{n}x_{i}} \right) = 2\left( 1 - \dfrac{1}{\bar{x}} \right) \\ \end{align*}

So I'd really appreciate it if anyone can see if the estimator I've found doesn't look weird, or give me any ideas on a way to verify this.

Answer 1

A quick and--usually--easy way to verify an estimator is to apply it to simulated data. I will describe this approach in a way that generalizes to any estimator in any situation.

Begin by coding your estimator. Here is an R implementation. Its input is a sample in an array x. It outputs the estimated parameter $\hat p(x).$

estimator <- function(x) 2 * (1 - 1 / mean(x))

You will need to generate random datasets. This function creates independent random Binomial$(m,p)$ variables truncated to exceed the value of $k:$

rbinom.trunc <- function(n, m, p, k = 0) qbinom(runif(n, pbinom(k, m, p), 1), m, p)

The simplest check is whether, on average, (1) the estimator's value is close to the parameter value for a range of parameter values and (2) that it gets closer as the sample size increases. That needs a double loop, implemented below using the outer function in R, which takes care of running the following simulate function for a specified sample size n and parameter value p:

simulate <- Vectorize(function(n, p, n.sim) {
  mean(apply(matrix(rbinom.trunc(n * n.sim, 2, p), n), 2, estimator))
}, c("n", "p"))

The third argument n.sim is the number of samples to generate. These are placed into a matrix, one sample per column, and the estimator is applied to each column to produce its estimate $\hat p.$ simulate returns the mean of all these estimates. (For a more detailed study, portray the entire set of estimates graphically with a histogram, probability plot, frequency plot, or whatever.) Here is an example of its use, where the values of $n,$ $p,$ and $\hat p$ are collected into a data frame for visualization:

n <- rev(c(2, 5, 10, 20, 50))
p <- c(0.05, 0.2, 0.5, 0.7, 0.9)
n.sim <- 5e2
X <- data.frame(n = factor(rep(n, length(p))), 
                p = rep(p, each = length(n)), 
                Estimate = c(outer(n, p, simulate, n.sim = n.sim)))

When I ran it, the output indicates the estimator is biased low for tiny values of $n,$ but once $n \ge 5$ or so, it is accurate. This doesn't mean it's a good estimator (we wouldn't expect the MLE to be a great one for small samples), but clearly it works.

This is a gussied-up version of the ggplot2 visualization of X created by ggplot(X, aes(p, Estimate, color = n)) + geom_point(). You can also make a usable, not-quite-so-pretty plot with the base plot command, as in with(X, plot(p, Estimate)).

---

Apart from the lines that specify the range of sample sizes n, the range of parameter values p, and the simulation size n.sim, this solution requires five lines of code for the five basic steps to estimate--generate--simulate--organize--summarize the results. Carrying out this kind of check often is so quick and easy (the computation time is negligible) that it's always worth doing when you care about your answer.

Answer 2

As a sanity check, you can think of this as an iid sample from a shifted Bernoulli distribution with parameter $q=p^2/(1-(1-p)^2)$. This gives you the MLE of $q$. You can then in turn use functional equivalence of MLEs to obtain the MLE of $p$.