Use law of total variance to find unconditional variance of overdispersed Poisson? First, I need to prove that the distribution of a RV X, where X|lambda ~ Pois(lambda), and lambda ~ gamma(a, B), is a negative binomial. I know that it is, but why negative binomial instead of another 2 parameter distribution? How do I prove negative binomial is the best alternative in the circumstance of an overdispersed Poisson, basically?
Then, I have to calculate V(X). I'm wondering if V(X) will help prove that negative binomial is best to use as an alternative in an overdispersed Poisson??
I looked at this post, and I desperately want to understand @probabilityislogic's math. I am new to notation associated with probability concepts (not that new, though), and I don't know how he/she jumps from line 2 to 3:

$λ_i∼Gamma(α,β)$
  On doing the integration/mixing over λi, you have:
  $Y_i(t_i)|αβ∼NegBin(α,p_i)$ where $p_i=\frac{t_i}{t_i+β}$

 A: 
How do I prove negative binomial is the best alternative in the circumstance of an overdispersed Poisson, basically?
I'm wondering if V(X) will help prove that negative binomial is best to use as an alternative in an overdispersed Poisson??

I really can't respond to those parts because it's unclear what you mean by "best" here. I don't think this kind of question is really answerable. If you have a gamma mixture of Poissons, the marginal distribution is negative binomial. If you don't, it may or may not be better to use it than overdispersed Poisson. It's often a good choice but that doesn't necessarily make it universally best over any unspecified circumstance and over every possible criterion of what you might want to be best at.
Responding to the last part (understanding probabilityislogic's mathematics) --
This I'll give an outline for. See also the Wikipedia page on the negative binomial distribution for a derivation in a slightly differently formulated but essentially identical problem.
Consider a case similar to the one done by probabilityislogic:
$$Y_i|\lambda_i\sim \text{Poisson}(\lambda_i)$$
$$\lambda_i\sim \text{Gamma}(\alpha,\beta)$$
Being slightly loose with notation for a moment:
$$f(y_i,\lambda_i)\,=p(Y_i=y_i|\lambda_i)\cdot f(\lambda_i)$$
To get the marginal for $Y_i$:
$$f_{Y_i}(y_i)=\int_0^\infty f(y_i,\lambda_i)d\lambda_i$$
$$=\int_0^\infty p(Y_i=y_i|\lambda_i)\cdot f(\lambda_i)\, d\lambda_i$$
So then just substitute the pmf and density inside, take constants (i.e. anything not in $\lambda_i$) out the front, collect like terms (leaving a $\lambda_i^\text{something}$ and an $\exp(-\lambda_i\cdot\text{something else})$ inside the integral).
Now we can play "spot the density" and recognize there's (apart from the required constant) a well known density inside the integral. Supply the constant (i.e. multiply by the required constant and divide by it out the front of the integral to keep the thing unchanged overall) and replace the integral by $1$. What's left should (once you see what terms need to be called "$p_i$") be a negative binomial for $Y_i$.
Beware: how the last part proceeds depends on whether you used the shape-rate or the shape-scale parameterization for the gamma density.
