Moment Estimation on Geometric Distribution I recently had trouble calculating the moment estimates for the parameter $p$ of the geometric distribution:
$$
P(X=k)=(1-p)^{k-1}p,\quad k=1,2,\cdots
$$
We know that there are two kinds of moment estimates for the geometric distribution parameter $p$. Generally we use a low-order form.
Assuming the sample mean is $\overline X$, it is not difficult to calculate a low-order moment estimate of the parameter $p$ as $\frac{1}{\overline X}$. i.e.$\hat{p}=\frac{1}{\overline X}$.
By calculation, the mathematical expectation of $\hat{p}$ is
$$
E(\hat{p})=np^n\sum_{k=n}^\infty \frac1k C_{k-1}^{n-1} (1-p)^{k-n},\quad \text{ where } 0<p<1
$$
In fact, $E(\hat{p})<\infty$, because the series on the right side of the equation can be proved to be convergent by D'Alembert's discriminant method.
Now I want to verify whether $\hat{p}$ is an unbiased estimator of $p$. Naturally, I need to calculate the limit of $E(\hat{p})$.
But through Mathematica's calculations, I found that the series on the right can be expressed as a hypergeometric2F1 function, and Mathematica cannot calculate this limit.
Sum[(k - 1)!/((n - 1)! (k - n)! k) (1 - p)^(k - n), {k, n, \[Infinity]}]

Hypergeometric2F1[n, n, 1 + n, 1 - p]/n

So can this limit be calculated? Can the unbiasedness of this estimator be verified? Thank you in Advance for your help.
 A: If $X \sim \operatorname{Geometric}(p)$ with $$\Pr[X = x] = (1-p)^{x-1} p, \quad x \in \{1, 2, \ldots\},$$ then for an IID sample $(X_1, \ldots, X_n)$, the sample total $$n \bar X = \sum_{i=1}^n X_i$$ has a negative binomial distribution $$\Pr[n \bar X = x] = \binom{x-1}{n-1} (1-p)^{x-n} p^n, \quad x \in \{n, n+1, \ldots \}.$$  The expectation of $1/(n \bar X)$ is the sum $$\operatorname{E}[1/(n \bar X)] = \sum_{x=n}^\infty \frac{1}{x} \binom{x-1}{n-1} (1-p)^{x-n} p^n,$$ which does not have an elementary closed form solution for general $n$.  So the expectation of the estimator $\hat p$, which is given by $$\operatorname{E}[\hat p] = n \operatorname{E}[1/(n \bar X)] = n p^n \sum_{x=0}^\infty \frac{1}{x+n} \binom{x+n-1}{n-1}(1-p)^x,$$ is expressed in Mathematica as $$\operatorname{E}[\hat p] = {}_2 F_1(1,1;1+n;1-p).$$  The fact that this is a function of $n$ should already be a clue that it does not equal $p$; thus $\hat p$ is biased.  If you just consider the small case $n = 2$, we get $$\operatorname{E}[\hat p \mid n = 2] = \frac{2 p (1-p+p \log p)}{(1-p)^2} \ne p.$$
