Friday, June 27, 2014

Probability density of sum of gamma and uniform

Probability density of sum of gamma and uniform

The probability density of the sum of gamma and uniform variables is the convolution of their densities, as it was for the sum of Gaussian and uniform variables, which we considered in a previous post. As before I'll define a simplifying unit boxcar function and a general boxcar to represent the uniform density.

matchdeclare (xx, numberp) $
tellsimp (U1(xx), if 0 <= xx and xx < 1 then 1 else 0) $
U(x, a, b) := 1/(b - a) * U1((x - a)/(b - a)) $

Let's let h be the density of a gamma variable, with parameters a and b.

h (x, a, b) := b^a / gamma(a) * x^(a - 1) * exp (-b*x);
h (x,a,b) := ( ba Γ (a) x a-1 exp ( ( -b ) x ) )

Here's the convolution integral. I'll quote the density functions to make the structure more obvious and therefore make further manipulation easier.

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

As in a previous blog post, I'll define some simplification rules in order to manipulate the integral to the point that Maxima's built-in integrate function can handle it. There's a slightly obscure strategem in the definition of the second rule. yy_minus_xx is assigned before turning off the simplification flag, since it is important that the rule match the simplified form of y - x, which is internally represented as y + (-1)*x.

matchdeclare (xx, symbolp, [aa, bb], all) $
matchdeclare (ee, all) $
matchdeclare ([ll1, ll2], all) $
yy_minus_xx : yy - xx $
simp : false $
tellsimp
('integrate (ee * 'U(xx, aa, bb), xx, ll1, ll2),
1/(bb - aa) * 'integrate (ee, xx, max (ll1, aa), min (ll2, bb))) $
tellsimp
('integrate ('U(''yy_minus_xx, aa, bb)*ee, xx, ll1, ll2),
'integrate ('U(xx, yy - bb, yy - aa)*ee, xx, ll1, ll2)) $
simp : true $

Now let's allow Maxima to apply the rules we just defined. The effect of quote-quote is just as if we just typed in the value of I.

I1 : ''I;
a 2 min (b 2 ,r) h ( r-s ,a 1 ,b 1 ) s b 2 -a 2

I originally defined the integral with h quoted in the integrand. Let's verbify it.

I2 : I, h;
b1a 1 a 2 min (b 2 ,r) ( r-s ) a 1 -1 - b 1 ( r-s ) s ( b 2 -a 2 ) Γ (a 1 )

I'll make a change of variables from r - s to u which simplifies the integrand somewhat.

I3 : changevar (I2, u - (r - s), u, s);
- b1a 1 r-a 2 r-min (b 2 ,r) u a 1 -1 - b 1 u u ( b 2 -a 2 ) Γ (a 1 )

Now we're ready to evaluate the integral. I'll assert some inequalities to forestall questions from 'integrate'. These are all necessary conditions for the definition of the gamma and uniform density functions.

assume (a[1] > 0, b[1] > 0, b[2] > a[2]) $
assume (r > a[2]) $
assume (min (b[2], r) > a[2]) $
I4 : I3, nouns;
- b1a 1 ( Γ (a 1 , b 1 r - b 1 a 2 ) b1a 1 - Γ (a 1 , b 1 r - b 1 min (b 2 ,r) ) b1a 1 ) ( b 2 -a 2 ) Γ (a 1 )

Terrific, now we have computed the convolution symbolically. Note the gamma function, Γ(a), and the upper incomplete gamma function, Γ(a, z), in the result. We'll take a look at a specific example, as given by specific values for the parameters.

I5 : subst ([a[1] = 1/4, b[1] = 1/2, a[2] = 5, b[2] = 12], I4);
- 2 1 4 Γ ( 1 4 , r 2 - 5 2 ) - 2 1 4 Γ ( 1 4 , r 2 - min (12,r) 2 ) 72 1 4 Γ ( 1 4 )

Let's define a function which has those parameter values baked in, and plot it.

conv_symbolic (r) := ''I5 $
plot2d (conv_symbolic, [r, 0, 20]) $
Produced by GNUPLOT 4.4 patchlevel 3 conv_symbolic r 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0 5 10 15 20

It's interesting that the convolution function is asymmetric. I guess that's not surprising, since the gamma density is asymmetric too. Well, that's all for now.

Thursday, June 26, 2014

Probability density of sum of Gaussian and uniform

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);
- g ( r-s ,a 1 ,b 1 ) U (s,a 2 ,b 2 ) s

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;
a 2 b 2 g ( r-s ,a 1 ,b 1 ) s 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;
πb 1 erf ( 2r - 2a 2 - 2a 1 2b 1 ) 2 - πb 1 erf ( 2r - 2b 2 - 2a 1 2b 1 ) 2 2πb 1 ( b 2 -a 2 )

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

I3 : ratsimp (I2);
- erf ( 2r - 2b 2 - 2a 1 2b 1 ) -erf ( 2r - 2a 2 - 2a 1 2b 1 ) 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 $
- erf ( 2r - 32 ) -erf ( 2r +2 5 2 ) 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]) $
Produced by GNUPLOT 4.4 patchlevel 3 conv_symbolic r 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 -10 -5 0 5 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])$
Produced by GNUPLOT 4.4 patchlevel 3 r discrete1 conv_symbolic -0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 -10 -5 0 5 10