# Computing the PDF of a quadratic function of two random variables

Given the function

$\mathcal{M} = g + Ah + Bh^2$

where $A$ and $B$ are constants and $g$ and $h$ are random variables with their distributions $f_G(g)$ and $f_H(h)$ known, is it possible to compute the probability density function $f_\mathcal{M}$?

We can obviously generate samples of it by creating realizations of $g$ and $h$ and then computing $\mathcal{M}$ (and then possibly form the empirical CDF, etc.), but I'm wondering if it is possible to "write down" an expression for $f_\mathcal{M}$.

A few notes:

1. If it is helpful, we can assume that $g$ and $h$ have normal distributions.
2. I'd be very interested in a solution that could accommodate the more complex case of $\mathcal{M} = g + \sum_i^n A_i h_i + B_i h_i^2$

Using the assumption of normality, we can get an expression for this in Maple, but it's not in closed form (it still has an integral sign in it). Here's what I did:

with(Statistics):
h := RandomVariable(Normal(mu[H], sigma[H]));
g := RandomVariable(Normal(mu[G], sigma[G]));
PDF(g + A * h + B * h^2, z) assuming A > 0, B > 0;


Without the assumptions on A and B I get an error, but with them, I get:

$$\int _{-1/4\,{\frac {{A}^{2}}{B}}}^{\infty }\!1/2\,\sqrt {2}{{\rm e}^{-1/2\,{\frac { \left( z-{\it \_t}-\mu_{{{\it G}}} \right) ^{2}}{{\sigma_{{{\it G}}}}^{2}}}}} \left( 1/2\,\sqrt {2}{{\rm e} ^{-1/8\,{\frac { \left( A+\sqrt {{A}^{2}+4\,B{\it \_t}}+2\,\mu_{{{\it H}}}B \right) ^{2}}{{B}^{2} {\sigma_{{{\it H}}}}^{2}}}}}{\frac {1}{\sqrt {\pi }}}{\sigma_{{{\it H}}}}^{-1}{\frac {1}{\sqrt {{A}^{2}+4\,B{\it \_t}}}}+1/2\,\sqrt {2}{{\rm e}^{-1/8\,{\frac { \left( A-\sqrt {{A}^{2}+4\,B{\it \_t}}+2\,\mu_{{{\it H}}}B \right) ^{2}}{{B}^{2}{\sigma_{{{\it H}}}}^{2}}}}}{\frac {1}{\sqrt { \pi }}}{\sigma_{{{\it H}}}}^{-1}{\frac {1}{\sqrt {{A}^{2}+4\,B{\it \_t}}}} \right) {\frac {1}{ \sqrt {\pi }}}{\frac {1}{\sqrt {{\sigma_{{{\it G}}}}^{2}}}}{d{\it \_t}}$$

(This is in an internal development version, so your results may be different. Full disclosure: I work for these guys.)

Actually, what might work even better is to see (by completing the square) that $$A h + B h^2 = B \left(\left(h + \frac{A}{2B}\right)^2 - \frac{A^2}{4 B^2}\right).$$ Suppose $h$ is normally distributed with parameters $\mu_H$ and $\sigma_H$, and let $C = \frac{A}{2B \sigma_H}$. Then so is $h /\sigma_H + C$ is normally distributed with variance 1 and mean $\mu_H/\sigma_H + C$. Its square is then distributed as a noncentral chi squared distribution with parameters $k = 1$ (a.k.a. $\nu = 1$) and $\delta = \mu_H/\sigma_H + C$. So we define $\tilde h = (h / \sigma_H + C)^2$, which is distributed in this way, and we examine $$A h + B h^2 = B \sigma_H^2((h+\frac{A}{2B})^2/\sigma_H^2 - \frac{A^2}{4B^2\sigma_H^2}) = B \sigma_H^2(\tilde h + C^2).$$ Now the PDF of that can easily be derived from the PDF of $\tilde h$. Unfortunately, that PDF involves a nasty hypergeometric function (a modified Bessel function of the first kind), so it doesn't make the evaluation of the PDF too much easier. However, it might make the theoretical analysis easier...

• The first part of your answer is similar to how we've attempted to solve this problem, though we had to solve it by hand (good to know Maple can solve this). I'll post what we came up with soon - it is a bit different. – Barron May 15 '12 at 17:16
• The second part is also interesting. @Emre suggested a convolution of the Gaussian part with the Noncentral Chi square, but I'm concerned with the lack of independence of the terms. Your separation of $g$ and $h$ might enable us to convolve $g$ with the PDF of $B\sigma_H^2(\tilde h + C^2)$. – Barron May 15 '12 at 17:22

First of all, we can decorrelate the linear and quadratic terms by completing the square:

$g + \sum A_i h_i + B_i h_i^2 \equiv g' + \sum B_i (h_i')^2$

The sum of squared nonstandard normal random variables appears to have no name or neat density, so I will compute it numerically, assuming that all the random variables are independent.

## Example

Remember that the density of the sum of two random variables is the convolution of their densities, and that convolution is a Fourier/Laplace transform pair with multiplication. Thus, we can easily convolve multiple pdfs by multiplying their Laplace- or Fourier transforms.

Since Pedro gave you an example in MATLAB, let me provide a complementary one in Mathematica. This toy example only has one quadratic term, but extension is easy.

mu = RandomReal[{0, 20}, 2]
sigma = RandomReal[{1, 5}, 2]
weights = RandomReal[{0, 1}, 2]
list1 = Table[PDF[NormalDistribution[mu[[1]], sigma[[1]]], x/weights[[1]]]/
weights[[1]], {x, 0, 50, 0.1}]
list2 = Table[
PDF[NoncentralChiSquareDistribution [1, mu[[2]]^2],
x/(sigma[[2]] weights[[2]])]/(sigma[[2]] weights[[2]]), {x, 0, 50, 0.1}];
ListPlot[{list1, list2}, PlotRange -> {{0, 300}, {0, 0.4}}]


This plot shows the discretizations of the normal and non-central chi-square random variables:

ListPlot[InverseFourier[ Times[##] &@
Apply[Sequence, Fourier[#] & /@ {list1, list2}]],
PlotRange -> {{0, 400}, {0, 0.06}}]


This plot shows the convolution; the final result:

• While $g$ and $h$ are independent, [$g + \sum A_i h_i$] and $Bh_i^2$ are not (because they both involve $h$). Would this be problematic for the convolution you propose? – Barron May 15 '12 at 13:51
• You are right; they need to be decorrelated. I made some corrections. – Emre May 15 '12 at 18:20

For the simple case of $\mathcal M = g + Ah$, the result is a convolution of the PDFs $f_G$ and $f_H$, i.e.

$$f_{\mathcal M}(x) = \int_{-\infty}^\infty f_G(t)\bar{f}_H(x-t)\,\mbox{d}t$$

where $\bar{f}_H$ is $f_H$ normalized such that $A$ disappears.

Multiple sums are convolutions of convolutions. This can quickly get quite messy.

One way of computing them numerically (in Matlab) is to use Chebfun (just as a disclaimer, I'm a member of the Chebfun development team):

>> normdist = @(m,s) @(x) 1/(s*sqrt(2*pi))*exp(-0.5*((x-m)/s).^2);
>> G = chebfun( normdist(0,1) , [-10,10] )
G =
chebfun column (1 smooth piece)
interval       length   endpoint values
[     -10,      10]       93  7.7e-23  7.7e-23
vertical scale = 0.4
>> H = chebfun( normdist(2,0.5) , [-10,10] )
H =
chebfun column (1 smooth piece)
interval       length   endpoint values
[     -10,      10]      169 6.7e-126  2.1e-56
vertical scale = 0.79
>> M = conv(G,H)
M =
chebfun column (2 smooth pieces)
interval       length   endpoint values
[     -20,       0]       59  1.2e-34    0.072
[       0,      20]       70    0.072  9.1e-35
Total length = 129   vertical scale = 1.6


We can verify this result by comparing it to the exact result

>> Mex = chebfun( normdist( 2 , sqrt(1^2+0.5^2) ) , [-20,20] )
Mex =
chebfun column (1 smooth piece)
interval       length   endpoint values
[     -20,      20]      154    3e-85  1.9e-57
vertical scale = 0.36
>> norm(M-Mex,inf)
ans =
5.5569e-16


I'm not using an infinite domain, since this seems to be broken in my current version of Chebfun. The result seems to be correct to about machine precision though.

• It isn't clear to me if/how we might incorporate the $B h^2$ term into this. Could that be added? – Barron May 15 '12 at 17:25
• @Barron: You could try generating the PDF of $h^2$ via the integral $H^2(x) \int_{-\infty}^\infty H(t)*H(x/t)\,\mbox{d}t$ and summing it as above. This is no longer a convolution, so in Chebfun you'd have to do something like H2 = chebfun( @(x) ... , [-10,10] ) where ... computes the above integral for some x. – Pedro May 15 '12 at 22:43