Repeated convolution of a continous uniform distribution



On Thu, 18 Jun 2009 11:30:03, weaker at directbox.com <weaker at directbox.com> wrote:

> I want to do a repeated convolution of a continous uniform distribution in a
> limited interval [-a,a].

A problem which is near and dear to my heart ...
I've appended a script I wrote to solve this problem
via pattern matching. Look for terms like y^n*U(y - a)*U(b - y)
and integrate them by hand.

For the 4-fold convolution I get:
(x^3*U(x)+(-4*x^3+12*x^2-12*x+4)*U(x-1)+(6*x^3-36*x^2+72*x-48)*U(x-2)
         +(-4*x^3+36*x^2-108*x+108)*U(x-3)+(x^3-12*x^2+48*x-64)*U(x-4))
 /6

8-fold is somewhat messier.
I think I went as far as 16-fold just for the fun of it one time.

Hope this helps. It was certainly entertaining to work on it.

best

Robert Dodier

PS. here's the script:

/* Copyright 2009 by Robert Dodier.
 * I release this work under terms of the GNU General Public License.
 */
matchdeclare (aa,freeof(y), bb, freeof(y), nn, freeof(y), cc, freeof(y));
defmatch (trunc_mono, cc*y^nn*U(y - aa)*U(bb - y), y);

hack (expr) := block ([m : length (expr), ss : 0],
    for k thru m do
        block ([term : part (expr, k), l],
            l : trunc_mono (term, x),
            new_term : cc*(bb^(nn + 1) - aa^(nn + 1))*U(bb - aa)/(nn + 1),
            /* print (k, l, new_term), */
            ss : ss + new_term),
    return(ss));

f1(x) := U(x) - U(x-1);

expr: expand (f1(x)*f1(y - x));

f2(y) := ''(hack (expr));

expr: expand (f2(x)*f2(y - x));

f4(y) := ''(hack (expr));

expr: expand (f4(x)*f4(y - x));

f8(y) := ''(hack (expr));

linel : 80;

grind (ratsimp (f2 (x)));
grind (ratsimp (f4 (x)));
grind (ratsimp (f8 (x)));

plot2d ([f2, f4, f8], [x, 0, 8], [plot_format, gnuplot]), U(x) := if x
>= 0 then 1 else 0;