# 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.

-
Have you considered that this mgf is the product of 2 mgfs? Given $X$ and $Y$ are independent random variables and $W = X+Y$, the mgf of $W$ is $M_W(t) = M_X(t) + M_Y(t)$. I am thinking something like $\frac{\lambda^2}{(\lambda-at)(\lambda-t)}=\begin{pmatrix}\frac{\lambda}{\lambda - at}\end{pmatrix}\begin{pmatrix}\frac{\lambda}{\lambda - t}\end{pmatrix}$ –  bryansis2010 Mar 27 '13 at 17:09
Try partial fractions, and then see if you can find what kind of random variable has an MGF of the form $\lambda/(\lambda-at)$ or $\lambda/(\lambda-t)$. –  Dilip Sarwate Mar 28 '13 at 3:28
@bryansis2010 If $X$ and $Y$ are not independent, what is the mgf of $W=X+Y$? –  Dilip Sarwate Mar 28 '13 at 3:30
@DilipSarwate One assumes that X and Y are independent here... :-) –  Did May 2 '13 at 23:59
@bryansis2010 don't you mean $M_w(t) = M_X(t)M_Y(t)$, i.e., the product? See Theorem 4.2.12 in Casella & Berger. –  chbrown Oct 29 '13 at 3:42

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


-
Why not using directly InverseLaplaceTransform? Is there some inner subtlety of the workings of Mathematica which would make preferable a detour by Fourier? –  Did May 3 '13 at 8:27
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.