Convolving random variables


#1

Hi all -

A while back I asked whether there is a package for convolving random variables in Julia (see original post). In the meantime, no package has been created and a few others have looked for similar functionality. As promised in the original post, I am sharing my code so that others may use it or perhaps improve upon it and create a proper package. (I also hacked together a way to represent the convolution as an object and with mathematical syntax (e.g. x = Normal(0,1) + Normal(0,1)). It might be an abuse of the language but I can share that if someone wants it.)

I am also hoping that someone who has a better understanding of signal processing and/or statistics might be able to show me how to can modify the function to convolve random variables for which the characteristic function is not available. My understanding is that there is a way to do this numerically with FFT. Unfortunately, this is outside my expertise and I have hit a wall. So any help would be much appreciated!

using Distributions,FFTW,Interpolations,Plots
function convolve(distributions,lb,ub,N = 2^9)
    i = range(0,length=N); dx = (ub-lb)/N
    x = lb .+ i*dx; dt = 2*π/(N*dx)
    c = -N/2*dt; t = c .+ i*dt; f = fill(1im+1.0,N)
    for dist in distributions
        @. f = f*cf(dist,t)
    end
    X = @. exp(-(0+1im)*i*dt*lb)*f
    fft!(X)
    v = @. dt/(2*π)*exp(-(0+1im)*c*x)*X
    #Use abs()to haldle very small negative values, e.g. -2.50e-12
    density = abs.(real(v))
    return interpolate((x,),density, Gridded(Linear()))
end
distributions = [Normal(0,1),Normal(0,1)]
lb,ub,=-6,6
convolved = convolve(distributions,lb,ub)
x = -5:.01:5
y = pdf.(Normal(0,sqrt(2)),x)
y′ = convolved.(x)
plot(x,y,color=:red,linewidth=2,grid=false,leg=false)
plot!(x,y′,color=:black,linestyle=:dash,linewidth=2)