Finding the nth moment of the geometric distribution: Why does interchanging the derivative and summation operators not work after n=1? I am trying (and failing) to find a recursive formula for the $nth$ moment of a geometric distribution. I have arrived at bogus results, and I think it has something to do with the convergence of the power series in question, and my own mathematically dubious operator-interchange. The problem:
Suppose we have a random variable $X$ that follows a geometric distribution (with parameter $p$, and let $q=1-p$), over the support $x \in \mathbb{N}$. The p.m.f is given by
$$
f(x) = pq^{x-1}
$$
I want to be able to find the $nth$ moment, that is, $E[X^n]$. I've seen the following clever trick used to find $E[X]$:
$$
E[X] = \sum_{x=1}^\infty xpq^{x-1} \\
= p\sum_{x=1}^\infty \frac{d}{dq}q^x \\
= p\frac{d}{dq}\sum_{x=1}^\infty q^x \\
= p\frac{d}{dq}(\frac{q}{1-q}) \\
= p\frac{(1-q) - (-q)}{p^2} = \frac{p}{p^2} = \frac{1}{p} = \frac{1}{1-q}
$$
This relies on interchanging the summation and differentiation operators (and happily factoring out the $p$, even though it is actually a function of $q$). I am trying to extend the idea to calculating the $nth$ moment (in terms of the $(n-1)$th moment):
First, notice that
$$
E[X^{n-1}] = \sum_{x=1}^\infty x^{n-1}pq^{x-1}
$$
And,
$$
E[X^{n}] = \sum_{x=1}^\infty x^{n}pq^{x-1}= \sum_{x=1}^\infty x^{n-1}pxq^{x-1} \\
= \sum_{x=1}^\infty x^{n-1}p\frac{d}{dq}q^{x}
$$
I imagine that the following step is the source of error:
$$
= \frac{d}{dq}\sum_{x=1}^\infty x^{n-1}pq^{x} \\
= \frac{d}{dq}q\sum_{x=1}^\infty x^{n-1}pq^{x-1} \\
= \frac{d}{dq}qE[X^{n-1}]
$$
This is obviously incorrect; setting $n=2$:
$$
E[X^2] = \frac{d}{dq}qE[X^{n-1}] = \frac{d}{dq}q\frac{1}{1-q} \\
= \frac{1}{(1-q)^2} = \frac{1}{p^2}
$$ and substituting into the formula for variance:
$$
\sigma^2 = E[X^2] - (E[X])^2
$$
yields:
$$
\sigma^2 = \frac{1}{p^2} - \frac{1}{p^2} = 0
$$
which I know is not true. Clearly I have made an error, and I don't know where, or more importantly, why. Can someone point me in the right direction please?

Starting form where you left off...
$$
E[X^{n}] = \sum_{x=1}^\infty x^{n}pq^{x-1}= \sum_{x=1}^\infty x^{n-1}pxq^{x-1} \\
= \sum_{x=1}^\infty x^{n-1}p\frac{d}{dq}[q^{x}]
$$
From this point, remember that p is a function of q, p = 1-q
$$
E[X^{n}] = \sum_{x=1}^\infty x^{n-1}(1-q)\frac{d}{dq}[q^{x}]
$$
Since (1-q) is a function of q and not X, we can factor this outside the sum but it must also remain outside the derivative operator. We can also interchange the summation and the derivative operator 
$$
E[X^{n}] = (1-q)\frac{d}{dq}[\sum_{x=1}^\infty x^{n-1}q^{x}]
$$
We then have,
$$
E[X^{n}] = (1-q)\frac{d}{dq}[q\sum_{x=1}^\infty x^{n-1}q^{x-1}] \\
= (1-q)\frac{d}{dq}[\frac{q}{(1-q)}\sum_{x=1}^\infty x^{n-1}(1-q)q^{x-1}] 
$$
Remember that 
$$
E[X^{n-1}] = \sum_{x=1}^\infty x^{n-1}(1-q)q^{x-1}
$$
This gives us
$$
E[X^{n}] = (1-q)\frac{d}{dq}[\frac{q}{(1-q)}E[X^{n-1}]]
$$
Now plugging in 2 for n we get, 
$$
E[X^2] = (1-q)\frac{d}{dq}[\frac{q}{(1-q)}E[X]] \\
= p\frac{d}{dq}[\frac{q}{(1-q)}\frac{1}{(1-q)}] \\
= p\frac{d}{dq}[\frac{q}{(1-q)^2}] \\
= p(\frac{1}{(1-q)^2} + \frac{2q}{(1-q)^3}) \\
= p(\frac{1 - q + 2q}{(1-q)^3}) = p(\frac{1 + q}{p^3}) \\ 
= \frac{1 + (1-p)}{p^2} = \frac{(2-p)}{p^2}
$$
Now substituting into the formula for variance:
$$
\sigma^2 = E[X^2] - (E[X])^2
$$
We get
$$
\sigma^2 = \frac{(2-p)}{p^2} - \frac{1}{p^2} = \frac{(1-p)}{p^2} \\
= \frac{q}{p^2}
$$
Which yields the correct variance formula.
I hope this helped :)
