Variance of Negative Binomial Distribution (without Moment Generating Series) Given the discrete probability distribution for the negative binomial distribution in the form 
$$P(X = r) = \sum_{n\geq r} {n-1\choose r-1} (1-p)^{n-r}p^r$$
It appears there are no derivations on the entire www of the variance formula $V(X) = \frac{r(1-p)}{p^2}$ that do not make use of the moment generating function. 
I have successfully managed to compute the mean without this as follows;
\begin{align*}
\mu = \sum_{n\geq r} n{n-1\choose r-1} (1-p)^{n-r}p^r 
&= \sum_{n\geq r} \frac{n(n-1)!}{(r-1)!(n-r)!}(1-p)^{n-r}p^r \\
&= \frac{r}{p} \sum_{n\geq r} \frac{n!}{r!(n-r)!}(1-p)^{n-r} p^{r+1}\\
\end{align*}
Having already factored our claimed mean of $r/p$, it remains to show that $\sum_{n\geq r} \frac{n!}{r!(n-r)!}(1-p)^{n-r} p^{r+1} = 1$ which is done by reindexing (both $r$ and $n$) and realizing this as the sum of a probability mass function for a negative binomial distribution. Indeed, letting $k = r+1$ followed by $m = n+1$, we find
\begin{align*}
\sum_{n\geq r} \frac{n!}{r!(n-r)!}(1-p)^{n-r} p^{r+1} &= \sum_{n\geq k-1}\frac{n!}{(k-1)!(n-k+1)!}(1-p)^{n-k+1}p^k\\
&=  \sum_{m\geq k}\frac{(m-1)!}{(k-1)!(m-k)!}(1-p)^{m-k}p^k\\
&=  \sum_{m\geq k}{m-1\choose k-1}(1-p)^{m-k}p^k  = 1
\end{align*}
Does anyone know of a way to demonstrate that $\sigma^2 = V(X) = \frac{r(1-p)}{p^2}$ in this fashion?
 A: No need to index twice. The process is quite similar to the way you get $E[x]$
I believe the problem here is how to get $E[x^2]$
Please refer to Markus Scheuer's efforts: 
In his answer he derived:
$$
\begin{align*}
E(x)&=rp^r\sum_{k=0}^{\infty}\binom{k+r}{k}(1-p)^k\\
\end{align*}.
$$
Follow the similar way, and apply to $E[X^2]$
$$
\begin{align*}
E(x^2)&=rp^r\sum_{k=0}^{\infty}(k+r)\binom{k+r}{k}(1-p)^k\\
\end{align*}
$$
We need a way to change $(k+r)$ to $(k+r-1)$, then the outside $(k+r)$ can be combined with $(k+r-1)!$ as $(k+r)!$, so I refer to the wiki page
$$
\begin{align*}
\binom{n}{k}=\binom{n-1}{k-1}+\binom{n-1}{k}\\
\end{align*}
$$
$$
\begin{align*}
(k+r)\binom{k+r}{k}&=(k+r)\binom{k+r-1}{k-1}+(k+r)\binom{k+r-1}{k}\\
&=(r+1)\binom{k+r}{k-1}+r\binom{k+r}{k}\\
\end{align*}
$$
So:
$$
\begin{align*}
E(x^2)&=rp^r\sum_{k=0}^{\infty}(k+r)\binom{k+r}{k}(1-p)^k\\
&=rp^r\sum_{k=0}^{\infty}[(r+1)\binom{k+r}{k-1}+r\binom{k+r}{k}](1-p)^k\\
&=rp^r\sum_{k=0}^{\infty}[(r+1)\binom{k+r}{k-1}](1-p)^k+rp^r\sum_{k=0}^{\infty}[r\binom{k+r}{k}](1-p)^k\\
\end{align*}
$$
then, still, use Markus Scheuer's idea in 
remember to use 
binomial series expansion
Finally, you will get
$$
\begin{align*}
E(x^2)&=\frac{r(1-r-p)}{p^2}\\
\end{align*}
$$
$$
\begin{align*}
V(X)&=E(x^2)-[E(x)]^2\\
&=\frac{r(1-r-p)}{p^2}+\frac{r^2}{p^2}\\
&=\frac{r(1-p)}{p^2}
\end{align*}
$$
A: This is too long for a comment, so I have it here as an answer.
Funny you ask this, since I was trying to figure this out yesterday. To prove that the Negative Binomial PDF does sum over $\mathbb{Z}_{\geq 0}$ to give $1$, you will need to make use of the binomial theorem for negative exponents (as Alex has indicated) and the fact posted at Negative binomial coefficient (but note the way this is written is for the "other" negative binomial distribution, with $K = X-r$).
The second moment $\mathbb{E}[X^2]$ is a bit tedious to compute. You will need to do reindexing twice rather than once with $\mathbb{E}[X]$. Unfortunately, the form of your negative binomial PDF is different from the one I worked with ($K = X-r$, as indicated above), so I don't have a sketch of this.
A: The negative binomial distribution with parameter $r$ is the distribution of  the number of times, $X$,  a Bernoulli experiment $B$ with probability $p$ has to be repeated independently to have it succeed for the $r$-th time.
Define $X_i$ to be the random variable denoting the number of times $B$ has to to be performed to succeed for the $i$-th time after having succeeded $i-1$ times. $X_i$ is a geometric random variable with probability of success $p$. Therefore the variance of $X_i$ is $\dfrac{1-p}{p^2}$. The $X_i$'s are all independent and hence we have,
$$ \operatorname{Var}[X]=\operatorname{Var} \left[ \sum_{i=1}^r X_i\right]=\sum _{i=1}^r \operatorname{Var}[X_i] = r\cdot\frac{1-p}{p^2}$$
A: I do like The Cryptic Cat's answer. I was also trying to find a proof which did not make use of moment generating functions but I couldn't find a proof on the internet. So I made an attempt.
Given
$$P(X = n) = \sum_{n\geq r} {n-1\choose r-1} (1-p)^{n-r}p^r,$$
where $n=\text{number of trials}$ and $r=\text{number of successes}$.
You have already to managed prove that
$$E(X) = \mu = \frac{r}{p}$$
Since expectation is a linear operator, rather than computing $E(X^2)$, it will be easier to find
$$\begin{align}
E(X(X+1)) &= E(X^2+X) = E(X^2)+E(X) \\
&= \sum_{n\geq r} n(n+1){n-1\choose r-1} (1-p)^{n-r}p^r \\
&= \sum_{n\geq r} \frac{n(n+1)(n-1)!}{(r-1)!(n-r)!}(1-p)^{n-r}p^r \\
&= \sum_{n\geq r} \frac{(n+1)!}{(r-1)!(n-r)!}\frac{r(r+1)}{r(r+1)}(1-p)^{n-r}p^r \\
&= \frac{r(r+1)}{p^2} \sum_{n\geq r} \frac{(n+1)!}{(r+1)!(n-r)!}(1-p)^{n-r}p^{r+2} \\
&= \frac{r(r+1)}{p^2} \sum_{n\geq r} {n+1\choose r+1}(1-p)^{n-r}p^{r+2}
\end{align}$$
Since
$\sum_{n\geq r} {n+1\choose r+1}(1-p)^{n-r}p^{r+2} = 1$,
$$E(X^2)+E(X) = \frac{r(r+1)}{p^2}$$
This implies that
$$E(X^2) = E(X^2)+E(X)-E(X) = \frac{r(r+1)}{p^2}-\frac{rp}{p^2} = \frac{r(1+r-p)}{p^2}$$
Therefore,
$$Var(X) = E(X^2) - [E(X)]^2 = \frac{r(1+r-p)}{p^2}-\frac{r^2}{p^2} = \frac{r(1-p)}{p^2}$$
A: $
\begin{split}
E(X^2)&=\sum \limits_{n=r}^\infty n^2\tbinom{n-1}{r-1}p^rq^{n-r}\\
&=\sum \limits_{t=0}^\infty (r+t)^2\tbinom{r+t-1}{r-1}p^rq^{t} \qquad (let \ n-r=t)\\
&=p^r\sum \limits_{t=0}^\infty (r+t)\frac{(r+t)!}{(r-1)!t!}q^{t}\\
&=p^r\sum \limits_{t=0}^\infty r\frac{(r+t)!}{(r-1)!t!}q^{t}+p^r\sum \limits_{t=0}^\infty t\frac{(r+t)!}{(r-1)!t!}q^{t}\\
&=r^2p^r\sum \limits_{t=0}^\infty \frac{(r+t)!}{r!t!}q^{t}+rp^r\sum \limits_{t=0}^\infty \frac{(r+t)!}{r!t!}tq^{t}\\
&=r^2p^rp^{-r-1} + rp^r\sum \limits_{t=0}^\infty \tbinom{-r-1}{t}t(-q)^{t}\\
&=\frac{r^2}{p}+rp^rq\sum \limits_{t=0}^\infty \tbinom{-r-1}{t}(-t)(-q)^{t-1}\\
&=\frac{r^2}{p}+rp^rq\sum \limits_{t=0}^\infty \tbinom{-r-1}{t}\frac{d(-q)^t}{dq}\\
&=\frac{r^2}{p}+rp^rq\frac{d}{dq}\sum \limits_{t=0}^\infty \tbinom{-r-1}{t}(-q)^t\\
&=\frac{r^2}{p}+rp^rq\frac{d(1-q)^{-r-1}}{dq}\\
&=\frac{r^2}{p}+\frac{(r+1)rq}{p^2}\\
Var(X)&=E(X^2)-[E(X)]^2\\
&=\frac{r^2}{p}+\frac{(r+1)rq}{p^2}-\frac{r^2}{p^2}\\
&=\frac{r(1-p)}{p^2}
\end{split}
$
