Calculate probability density function from moment generating function I have moment generating function 
$$ M_z (t) = \dfrac {\lambda^2} {(\lambda-at) (\lambda-t)}, $$ and I'm trying to calculate the PDF from this function. I feel like some kind of inverse transform is required here, but I can't figure it out.
 A: Inverting generating functions is not always easy. It does help to state any constraints on the parameters (which are left unstated here by the OP). Anyway, the approach is: First, convert your mgf into a characteristic function (i.e. replace t -> I t):

Next, invert the characteristic function to yield the pdf (using Mathematica here), using an inverse fourier transform. To get a 'neat' result, it appears we need the parameter constraints {a > 0, $\lambda$ > 0}, so I have added these as assumptions:

HeavisideTheta[x] is a function that is 0 for x<0, and 1 for x > 0, so our pdf is:
 
Finally, a quick check to make sure all is OK:  generate the mgf from the pdf to see it is the same as that we started with (here using the Expect function in mathStatica/Mathematica):

... and all is OK :)
Finally, here is a quick plot of the derived pdf, as parameter $\lambda$ varies:
PlotDensity[f /. {a -> 2, \[Lambda] -> {2,3,4}}]


A: (This explains in detail a suggestion by @Dilip then expands briefly on another one by @bryanis2010, both leading to excellent general approaches. Thus, the aim of this post is to provide a mathematical answer to a question which has none, so far.)
As every rational function, $M_Z$ can be decomposed into partial fractions, here, for every $a\ne1$,
$$
M_Z(t)=\frac{u}{\lambda-at}+\frac{v}{\lambda-t},
$$
for some parameters $u$ and $v$ one can readily identify as $u=-a\lambda/(1-a)$ and $v=\lambda/(1-a)$.
Now, one of the easiest to remember generating functions is that, if $X$ is standard exponential, then
$$
M_X(t)=E[\mathrm e^{tX}]=\int_0^{+\infty}\mathrm e^{tx}\mathrm e^{-x}\mathrm dx=\frac1{1-t}.
$$
Hence, for every positive $\mu$,
$$
\frac1{1-\mu t}=E[\mathrm e^{t\mu X}].
$$
Applying this to $\mu=1/\lambda$ and to $\mu=a/\lambda$ yields
$$
M_Z(t)=\frac{u}{\lambda}E[\mathrm e^{taX/\lambda}]+\frac{v}{\lambda}E[\mathrm e^{tX/\lambda}].
$$
This indicates that 
$$
M_Z(t)=\frac{u}{\lambda}\int_0^{+\infty}\mathrm e^{tax/\lambda}\mathrm e^{-x}\mathrm dx+\frac{v}{\lambda}\int_0^{+\infty}\mathrm e^{tx/\lambda}\mathrm e^{-x}\mathrm dx,
$$
that is,
$$
\int_0^{+\infty}\mathrm e^{tx}f_Z(x)\mathrm dx=-\frac{a}{1-a}\int_0^{+\infty}\mathrm e^{tx}\mathrm e^{-x\lambda/a}\frac{\lambda}a\mathrm dx+\frac1{1-a}\int_0^{+\infty}\mathrm e^{tx}\mathrm e^{-x\lambda}\lambda\mathrm dx,
$$
hence, finally,
$$
f_Z(x)=\frac{\lambda}{1-a}(\mathrm e^{-x\lambda}-\mathrm e^{-x\lambda/a}).
$$
The simplest method to solve the case $a=1$ might be to consider the limit when $a\to1$ of this expression, yielding the correct gamma density
$$
f_Z^{(a=1)}(x)=\lambda^2x\mathrm e^{-x\lambda}.
$$

Another method, suggested by @bryanis2010 in comments is to note that $M_Z$ is the product of two generating functions, each easily identifiable, namely
$$
M_Z(t)=\frac{\lambda}{\lambda-at}\cdot\frac{\lambda}{\lambda-t}=E[\mathrm e^{atX/\lambda}]\cdot E[\mathrm e^{tX/\lambda}],
$$
hence $Z$ is distributed as $(aX+Y)/\lambda$ where $Y$ is also standard exponential, independent from $X$. This yields $f_Z$ as the convolution of two well known densities, for the same end result as above, the $a=1$ case (gamma random variable) being particularly direct.
