5
$\begingroup$

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.

$\endgroup$
6
  • 5
    $\begingroup$ Where specifically are you looking for help? We're not in the homework-checking business, but we are glad to respond to specific questions. $\endgroup$
    – whuber
    Commented Oct 3, 2022 at 17:23
  • $\begingroup$ of course i know you are are not in the homework-checking business, you don't have to say such obvious things @whuber $\endgroup$
    – aliocha
    Commented Oct 3, 2022 at 19:31
  • 2
    $\begingroup$ Given that's exactly what you are asking for, the obvious appeared to be necessary. $\endgroup$
    – whuber
    Commented Oct 3, 2022 at 19:34
  • 1
    $\begingroup$ What I might like or don't like is of no importance. The site guidelines for what kinds of questions to ask and what is on topic are what matters: please see stats.stackexchange.com/help/dont-ask and stats.stackexchange.com/tags/self-study/info. $\endgroup$
    – whuber
    Commented Oct 3, 2022 at 20:19
  • 2
    $\begingroup$ +1. The updated question, asking for methods to verify a putative estimator, is excellent. Have you considered applying it to simulated data? $\endgroup$
    – whuber
    Commented Oct 3, 2022 at 21:47

2 Answers 2

6
$\begingroup$

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.

enter image description here

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.

$\endgroup$
1
  • 2
    $\begingroup$ what an excellent way to check the estimator, I hadn't thought of using simulation, thank you! $\endgroup$
    – aliocha
    Commented Oct 4, 2022 at 0:53
5
$\begingroup$

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$.

$\endgroup$
3
  • $\begingroup$ i was thinking about your suggestion, and you mean a new random variable $Y$, with this pmf: \begin{array}{lclll} \mathbf{P}(Y = y) & = & \left[ \dfrac{p^2}{(1-(1-p)^2)} \right]^{y} \cdot \left[ \dfrac{2p(1-p)}{(1-(1-p)^2)} \right]^{1-y} \cdot \mathbf{1}_{\{0,1\}}(y) \end{array} $\endgroup$
    – aliocha
    Commented Oct 4, 2022 at 23:11
  • 1
    $\begingroup$ Yes, if you define $y_i=x_i-1$, then the $y_i$'s are iid Bernoulli. Hence, $\hat q=\bar y=\bar x-1$. Expressing $p$ in terms of $q$ you'll find (after some algebra) that $p=2q/(1+q)$. Functional equivalence implies that the MLE of $p$ is the same function of the MLE of $q$, that is, $\hat p=2\hat q/(1+\hat q)$ which simplifies to your expression when you substitute $\bar x-1$ for $\hat q$. $\endgroup$ Commented Oct 5, 2022 at 13:29
  • $\begingroup$ Perfect. Thank you! $\endgroup$
    – aliocha
    Commented Oct 5, 2022 at 14:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.