Computing the PDF of a quadratic function of two random variables

Question

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$

Answer 1

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

Answer 2

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:

Answer 3

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.