I need help with some maths/programming problem...
I try to smoothe a discrete probability function (of a data set) with a Gaussian distribution. I understand that I have to use a convolution:
(g*f)(m) = sum(n=1 -> m) f(n)g(m-n); this is a sum in my case because it's a discrete convolution. The function f(x) is my original probability distribution, the function g(x) is a Gaussian distribution.
I'm getting stuck with how to actually make the convolution work... it seems easy, but I seem to be getting too big numbers.
More detail on my original distribution: it is defined for 0 <= x <= 6 in steps of 0.05 (I've got a concrete number at each step) and it integrates to an area of 82 (that's my total number of data objects). Now I want to convolve it with a Gaussian with sigma = 0.4, and I use
g(i) = 1/(0.4 * sqrt(2.0*3.1415926))*exp(-((x(i)**2)/(2* (0.4**2))))
which should be normalised to 1.
First question now: where does the centre of the peak of the Gaussian have to be - I think it shouldn't matter, should it; but my program spits out different numbers if I shift the peak. I think however that if I leave the peak at zero, then I only define half of the function, since the range I'm looking at goes from 0 to 6?
Even if I do this though (meaning: I only use half of a normalised Gaussian, which then has an area of 0.5), the area of my final convolved function is something like 857, way too high. What could be happening here?
For the actual convolution my code looks as follows; n_z = 121 ( because xmax = 6 in steps of 0.05), g(i) is the above definition, f(j) is my probability distribution, and int in the code is my attempt to calculate the area of the final result.
I'm also not sure about the i-j not being able to be below zero; that was just my thoughts because I only defined both functions for values greater than zero.
c Perform the convolution
do i = 1, n_z
do j = 1, i
c Avoid range violation: i-j cannot be below zero
if (i-j.le.0) then
conv(i) = conv(i) + (f(j) * 0.0)
conv(i) = conv(i) + (f(j) * g(i-j))
int = int + (conv(i)*0.05)
Any help appreciated a lot - thank you!!