MLE of Geometric distribution - consistency and variance of inverse arithmetic sum I want to calculate the MLE and its consistency of Geometric distribution:
$$\mathbb{P}(X=x)=p(1-p)^{x-1}$$
$$
l(p) = p(1-p)^{x_1-1}\cdot p(1-p)^{x_2-1} \cdot \cdots \cdot p(1-p)^{x_n-1} \\
L(p)=n\ln p + (x_1 + \cdots +x_n -n)\ln(1-p) \\
(L(p))'= \frac{n}{p} - \frac{x_1 + \cdots + x_n - n}{1-p} \\
\hat{p} = \frac{n}{x_1 + \cdots + x_n} = \frac{1}{\overline{X}}
$$
Now, to check the consistency I would use Markov's Inequality:
$$
\mathbb{P}(|\hat{p_n} - \mathbb{E}(\hat{p_n})| \geq \epsilon) \leq \frac{\operatorname{Var}(\hat{p_n})}{\epsilon^2}
$$
I am stuck on calculating Variance though. Variance of arithmetic sum would be a breeze - here though it is an inverse of it. How can i proceed? Is it just the inverse of $\operatorname{Var}(\overline{X})$ = $\frac{\operatorname{Var}(X)}{n}$ which would equal $\frac{n}{\operatorname{Var}(X)}$? What is the $\mathbb{E}(\hat{p_n})$?
 A: It is a general fact that maximum likelihood estimators are consistent under some regularity conditions. In particular these conditions hold here because the distribution of $X$ is a member of a regular exponential family. Some discussion can be found in these lecture notes, but you can also find these results on a textbook covering asymptotic theory.
Consistency can be justified from the law of large numbers, which says 
$$\overline X \stackrel{P}\longrightarrow \operatorname E\,[X_1]=\frac1p$$
By continuous mapping theorem, this implies
$$\frac1{\overline X} \stackrel{P}\longrightarrow p$$
Alternatively, you can use a Taylor expansion on moments to say that $\operatorname E\left[\frac1{\overline X}\right]\approx p$ and $\operatorname{Var}\left[\frac1{\overline X}\right]\approx \frac{p^2(1-p)}{n}$ for large $n$, so that $\operatorname E\left[\frac1{\overline X}\right]\to p$ and  $\operatorname{Var}\left[\frac1{\overline X}\right]\to 0$ as $n\to \infty$. This is a sufficient condition for convergence in probability.
A: Here is how one may compute the expectation of $\hat{p}_n$. Let $S_n:=\sum_{i=1}^n X_i$. Then $S_n-n\sim \text{NB}(n,1-p)$ (negative binomial distribution) and the expectation of $\hat{p}_n$ takes the following form:
\begin{align}
\mathsf{E}\hat{p}_n&= \sum_{k\ge 0}\frac{n}{k+n}\binom{k+n-1}{k}(1-p)^k p^n \\
&= p\cdot{}_2F_1(1,1;n+1;1-p),
\end{align}
where ${}_2F_1$ is the hypergeometric function. Since
$$
{}_2F_1(1,1;n+1;1-p)=n\int_0^1 \frac{(1-x)^{n-1}}{1-x(1-p)}\,dx,
$$
it follows that for $p\in (0,1)$,
$$
p< \mathsf{E}\hat{p}_n< \frac{p}{1-n^{-1}}.
$$
That is $\hat{p}_n$ is asymptotically unbiased, i.e., $\mathsf{E}\hat{p}_n\to p$ as $n\to\infty$.
A: MLE of $p$ for Geometric Distribution is biased. You can use Jensen's Inequality to show this:
$$\
E(\frac{1}{\overline{X}}) > \frac{1}{E(\overline{X})} > \frac{np}{n} > p
$$
Note that $\frac{1}{x}$ is a strictly convex function for $x \in (0, c)$ where $c$ is a constant which is bounded wagely, $c < \infty$. Hence, we can have strong inequality ($<$ instead of $\leq$). To compute the exact value you can use the fact that sum of n i.i.d Geometric r.vs follows Negative Binomial Distribution. Hence,
$$\
E(\frac{1}{\overline{X}}) = nE(\frac{1}{Y}) \text{ where, } Y \sim NegBin
$$
For consistency you can proceed like this:
$$\
P(|\frac{1}{\overline{X}}-p|\leq\epsilon)=P(\frac{1}{\epsilon +p} \leq \overline{X} \leq \frac{1}{p-\epsilon}) = P(\frac{n}{\epsilon +p} \leq Y \leq \frac{n}{p-\epsilon})\\ \text{ where, } Y \sim NegBin
$$
You can compute the above and check consistency. I guess it will turn out to be consistent.
Edit: It turns out the cdf of Negative Binomial Distribution is not straight forward to solve. We need to use some inequality.
A: A brief simulation in R will illustrate some of the issues
you need to consider. R uses a version of the geometric
distribution taking values $0, 1, 2, \dots.$ Your version
takes values $1, 2, 3, \dots,$ hence the $+1$ in the R code
to illustrate your version.
Let $n=10, p = 0.3:$
set.seed(329);  n = 10;  p = .3
a = replicate(10^6, mean(rgeom(n, p)+1))  # averages of 10 obs
mean(a);  var(a)        
[1] 3.334036              # aprx 1/p = 1/.3 
[1] 0.7815287             # aprx (.7/.3^2)/10 = .7778

The MLE is $\hat p_n = 1/\bar X,$ as you say. However, it is clearly not unbiased. (With a million iterations, we should get about
two place accuracy.)
p.est = 1/a
mean(p.est);  var(p.est)
[1] 0.3217807             # NOT aprx 0.3
[1] 0.007767821

Estimation of $p$ is an important problem for some medical applications, and it is not trivial. Wikipedia moves quickly to a discussion of Bayesian approaches. You can google to access
a few scholarly papers on the estimation of geometric $p$ that are
not behind pay walls.
