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.

No comments:

Post a Comment