# Sample from a sum of two RVs by convolution of their samples via fourier transform in Excel

I apologies if this is off topic for this site, but I am trying to figure out (and not having much luck) how to get an empirical distribution function of a sum of two random variables given two samples of equal length (say 256) from their respective distributions. I understand that the distribution is given by a convolutions of their distributions, which can be computed as a product of their characteristic functions.

So, to test, I take two time series of length 256 each from standard normal distribution, compute their Fourier transforms using Fourier transform tool from Data Analysis toolpack in Excel, multiply the results using complex multiplication function IMPRODUCT, and then apply inverse Fourier transform. However, the resulting time series does not have a variance of 2 (it is more like ~240) which should be the case for the sum of two standard normal RVs. It might be some normalization issue, but I am not sure how and why it comes up.

Also, I am not sure if I am supposed to be doing zero-padding of the two original time series (say +256 zeros for each). I have read that if one wants to perform regular convolutions (instead of a "circular" one, whatever that is), one have to do this padding with zeros. This padding, however, seems to make the final result also double the length, but I would like to get the estimated distribution of the sum of the same length as the inputs, so not sure what needs to be done here.

For example, given two identical samples of size 9 from uniform distribution

$$(0,0.125,0.250,0.375,0.50,0.625,0.750,0.875,1.0)$$,

their convolution would produce 81 pairs with values

$$(0,0.125,0.250,0.375,0.50,0.625,0.750,0.875,1.0,1.125,1.250,1.375,1.50,1.625,1.750,1.875,2.0)$$

with respective counts

$$(1,2,3,4,5,6,7,8,9,8,7,6,5,4,3,2,1)$$

which gives the triangular distribution. What I am looking to do, is to produce this without doing these convolution sums, but by using Fourier transform.

This is what Wikipedia says on the calculation of empirical characteristic function.

If $$X_{i},i=1,2,\ldots ,n$$ are i.i.d. observations, then the empirical characteristic function $$\varphi _{n}(t)$$ is defined as

$$\varphi _{n}(t)={\frac {1}{n}}\sum _{j=1}^{n}\exp(itX_{j}).\tag{1}$$

And this is what it says about application of the characteristic function to sums of random variables.

In particular, $$\varphi _{X+Y}(t) = \varphi _{X}(t) \varphi _{X}(t)$$, i.e. the characteristic function of the sum of RVs is the product of their respective characteristic functions.

So, from the above two, it seems that one can construct characteristic functions directly from the samples and multiply them together to get a characteristic function of the sums. The only bit left is to get back into the "sample space". That was my understanding.

A discreet Fourier transform (DFT) of a sequence is

$$\hat f_k=\sum _{j=0}^{n-1} f_j \exp(i2\pi k j/n).\tag{2}$$

So, if we take the probability of observing each value in the sample on $$n$$ observations the same, ie $$f_j = 1/n$$, then the difference between empirical characteristic function and DFT is that $$X_j$$ in the former is replaced with $$2\pi j /n$$ in the latter, so the characteristic function looks like some kind of "non-uniform" DFT. And this is where I seem to get stuck, I don't know how to compute this non-uniform DFT. Also, the discreet $$t$$ points of the transform need to be the same for both samples as I would need to take the product of the two empirical characteristic functions at those points.

• wait, you draw 256 samples from a distribution, and then act as if that is the distribution? It's not! It's a sample vector!! If you approximate the distribution by finding the histogram, you might be doing something sensible... (also, really, Excel is not the tool of choice. Try Python, R or Octave instead. Free tools that actually do what you need). Jul 21, 2017 at 14:34
• Thank you for your comment. If I am not mistaken, one can take each observation as coming from a discrete distribution with probability of 1/n (n being 256 in this case), so effectively sample is like a distribution. This is also what I saw being done for calculation of empirical characteristic function, just applied to the sample. (Not sure, however, if one needs to sort the sample first). Jul 21, 2017 at 14:36
• no, an observation is an observation, and not an estimate of the density.You need to calculate the empircal distribution first, and then you can convolve that. Jul 21, 2017 at 14:38
• As @MarcusMüller says, you're completely off-base here. Please read this question and its answers and revise your question so that it makes a bit more sense.
– Peter K.
Jul 21, 2017 at 14:46
• Your addition has ordered the samples... and they "happen" to be completely uniformly spaced. That won't happen with a truly random variable. The convolution of a 9-sample vector with itself will yield a 17-sample vector (well, non-zero if all 9 samples are non-zero), not 81. Again, you seem to be mixing the samples from the distribution with the underlying distribution.
– Peter K.
Jul 21, 2017 at 16:15

I still think you're completely confused about what you're trying to do, so I'm going to go back and assume your question is:

Given two data sets $x_n, n=0,\ldots, N-1$ and $y_n, n=0,\ldots, N-1$, how do I compute the probability density function (PDF) of $x_n + y_n$ using Excel?

Finding the PDF of $x_n$

Finding the empirical PDF of a data sequence is a little tricky. The most straightforward way is to find the histogram of the data and normalize it so that it sums to 1.

In pseudocode:

h_x = histogram(x)
pdf_x = h_x/sum(h_x)


It's not Excel, but this Google Sheet shows how that can be achieved.

The columns are:

• x : the random sequence.
• bin start : the start of each bin (next one is the end of this bin)
• count : the number of samples of x that fall between the bin start and bin end.
• pdf_hat : a stab at the empirical PDF using the histogram.

The same approach can be used for finding the PDF of $y$.

Finding the PDF of $x_n + y_n$

As you say, the PDF of the sum is the convolution of the two individual PDFs (provided $x_n$ and $y_n$ are independent). See section 7.2 here.

Unfortunately, Google Sheets doesn't have an FFT function, but you could do it element by element, which is what I've tried to show in the columns:

• index : index of the PDF estimates.
• Realization 1 : the pdf hat values for one realization.
• Realization 2 : the pdf hat values for another realization.
• Convolution : the element-by-element convolution.

The resulting graph is an estimate of the PDF of the sum, which approximates the expected triangle function quite well. How to do better?

Generally, a better way to estimate the PDF of a random variable is to use kernel density estimators. Once you've got the individual PDF estimates then you can just convolve them. Or, if you have the summed data, just apply the kernel density estimator directly to that data.

Also, as you say, a better way to do convolution may be to use the FFT. To avoid "time aliasing" you need to zero-pad your signals (PDFs in this case) to be at least $N+M-1$ long (where one signal is length $N$ and the other $M$).

• Thank you for taking time to read and reply to my post. I believe that the problem is "reduced" to finding a non-uniform DFT of a constant function $1/n$ sampled at the "times" given by the observed random numbers in the sample. This should results in an empirical characteristic functions of the generating distribution. Once two such empirical characteristic function are obtained, they can be multiplied together to get a characteristic function for a sum of two RVs, which can be then inverted to get empirical distribution of the sum. Jul 22, 2017 at 11:47
• Thank you for your answer. However, I am not looking to find the PDF of the sum - I want to obtain a sample. That is, given a sample $x_n$ and $y_n$ with $n = 0, ..., N-1$, I want to obtain a new sample $z_n$ of length $N$ which would correspond to a sample from $x_n + y_n$. There is a brute-force approach to finding the "full" sample of size $N^2$ by summing each element of $x_n$ with each element of $y_n$, but this is too computationally expensive for my purposes, so I am looking for a way to get a "sub-sample" and I thought that maybe FFT can be used for that. Mar 30, 2022 at 13:14