Probability density of sum of Gaussian and uniform

The probability density function of the sum of a Gaussian variable and a uniform variable is the convolution of their densities. Let's compute that convolution symbolically, and compare it with a numerical approximation.

Constructing the symbolic convolution

First let's define the Gaussian density g and the uniform density U. These depend on some parameters.

g(x, a, b) := 1/(b*sqrt(2*%pi))*exp(-(1/2)*(x - a)^2/b^2) $

U(x, a, b) := 1/(b - a) * U1((x - a)/(b - a)) $

I've defined the general uniform density in terms of a unit uniform density U1 defined on the interval [0, 1). (Maxima doesn't have a built-in notion of a unit step or boxcar function so I'll do this part myself.) I'll make U1 a simplifying function in order to prevent it from being evaluated except when the argument is a number.

matchdeclare (xx, numberp) $

tellsimp (U1(xx), if 0 <= xx and xx < 1 then 1 else 0) $

Now let's construct the convolution integral. In general the convolution conv(f, g) = integral f(r - s) g(s) ds. I'll quote g and U so that they are not evaluated; that makes it easier to work with the integral in this case.

I : integrate ('g(r - s, a[1], b[1]) * 'U(s, a[2], b[2]), s, minf, inf);

⇒

${\int }_{-\infty }^{\infty }g(r-s,a_1,b_1)U(s,a_2,b_2)ds$

Simplifying and evaluating the convolution integral

The convolution integrand contains a function which is nonzero only on an interval, and constant there, so we can pull out the constant factor and change the limits of integration. I'll define a simplification rule for that.

matchdeclare (xx, symbolp, [aa, bb], all) $

matchdeclare (ee, all) $

matchdeclare ([ll1, ll2], all) $

simp : false $

tellsimp ('integrate (ee * 'U(xx, aa, bb), xx, ll1, ll2), 1/(bb - aa) * 'integrate (ee, xx, max (ll1, aa), min (ll2, bb))) $

simp : true $

First step toward evaluation is to let Maxima apply the simplification rule I just defined.

I1 : ''I;

⇒

$\frac{{\int }_{a_2}^{b_2}g(r-s,a_1,b_1)ds}{b_2-a_2}$

Terrific, now it's just a Gaussian density in the integrand. It's a noun (i.e., quoted function), so verbify it. Then 'integrate' can compute its integral.

I2 : I1, nouns;

⇒

$\frac{\frac{\sqrt{\pi }{b}_1\mathrm{erf}\left(\frac{\sqrt{2}r-\sqrt{2}{a}_2-\sqrt{2}{a}_1}{2b_1}\right)}{\sqrt{2}}-\frac{\sqrt{\pi }{b}_1\mathrm{erf}\left(\frac{\sqrt{2}r-\sqrt{2}{b}_2-\sqrt{2}{a}_1}{2b_1}\right)}{\sqrt{2}}}{\sqrt{2}\sqrt{\pi }{b}_1\left(b_2-a_2\right)}$

Rearrange it somewhat, just to make it nicer to look at.

I3 : ratsimp (I2);

⇒

$-\frac{\mathrm{erf}\left(\frac{\sqrt{2}r-\sqrt{2}{b}_2-\sqrt{2}{a}_1}{2b_1}\right)-\mathrm{erf}\left(\frac{\sqrt{2}r-\sqrt{2}{a}_2-\sqrt{2}{a}_1}{2b_1}\right)}{2b_2-2a_2}$

I'll pick out some typical values for the parameters, and bake them into a function to evaluate the convolution.

I4 : subst ([a[1] = -6, b[1] = 1/2, a[2] = 2 , b[2] = 9], I3);

conv_symbolic (r) := ''I4 $

⇒

$-\frac{\mathrm{erf}\left(\sqrt{2}r-3\sqrt{2}\right)-\mathrm{erf}\left(\sqrt{2}r+2^{\frac{5}{2}}\right)}{14}$

Let's make a plot. The convolution of a Gaussian and a uniform looks like a softened boxcar. Each end of the boxcar is a shifted and rescaled erf.

plot2d (conv_symbolic, [r, -10, 10]) $

Numerical approximation to the convolution via FFT

Here I'll just paste in the discrete convolution via FFT which was the subject of a recent blog post.

load (fft) $

conv (a, b) :=

block

([n : length(a) + length(b) - 1],

inverse_fft (fft (zero_extend (a, next2(n))) * fft (zero_extend (b, next2(n)))),

realpart (next2(n) * %%),

rest (%%, n - next2(n))) $

zero_extend (a, n) := makelist (if i <= length(a) then a[i] else 0, i, 1, n) $

next2 (n) := block ([i : ilog2 (n)], if n = 2^i then n else 2^(i + 1)) $

ilog2 (n) := block ([i : 0], while n > 1 do (n : floor (n/2), i : i + 1), i) $

Comparing the symbolic and numerical results

I'll create a grid of values and evaluate the Gaussian and uniform densities on the grid. The parameters for the densities are the same as what were baked into conv_symbolic.

n : 100 $

dx : (10 - (-10))/n $

xx : makelist (-10 + (i - 1) * dx, i, 1, n), numer $

gx : dx * makelist (g (xi, a[1], b[1]), xi, xx), a[1] = -6, b[1] = 1/2, numer $

Ux : dx * makelist (U (xi, a[2], b[2]), xi, xx), a[2] = 2, b[2] = 9, numer $

Then I'll compute the convolution of those discretized densities.

conv_numerical : conv (gx, Ux) $

conv_numerical_density : conv_numerical / dx $

conv_support : makelist (-20 + (i - 1) * dx, i, 1, 2*n - 1) $

Finally I'll plot both the discretized approximation and the symbolic density. The agreement is pretty close. That's comforting.

plot2d ([[discrete, conv_support, conv_numerical_density], conv_symbolic], [r, -10, 10])$