Variance of Negative Binomial Distribution (without Moment Generating Series)

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*} $$

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

$ \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} $

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.

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