Recursive convolution of a PDF with itself

I was playing around with some simulations and I am trying to understand and produce an analytical solution for the behaviour I see.

I generated $$512$$ sets of data, each set containing $$10000$$ data points. All data sets were drawn from a Rayleigh distribution with the same shape parameter $$\sigma$$. I label each set $$X_{i}$$, which has elements $$X_{i} = \{x_{j},...,x_{M}\}$$, where $$M = 10000$$.

I then take the average of the data sets as $$\bar{X}_{N} = \frac{1}{N} \sum_{i}^{N} X_{i}$$ to be explicit, not the average of the $$X_{i}$$, but the average of the sets. The resulting average still contains $$M = 10000$$ points. Plotting the averaged data sets gives the histograms, as follows:

What we can clearly see is that as more and more data sets are included in the average, the more normally distributed the resultant data set becomes.

This is to be expected, I think, because of the central limit theorem. Now to my question:

I would like to show this mathematically. I know that the sum of random variables is described by the convolution from the parent distribution of those variables.

I imagine the result to be a two parameter distribution, $$F(\sigma, N)$$, $$\sigma$$ being the same shape parameter of the Rayleigh distribution. And $$N$$ the number of averages, or convolution iterations as per $$F(x)_{N} = \int_{-\infty}^{+\infty} F_{N - 1}(X) f(x - X) \ dX$$ where $$f(x)$$ is the original distribution of interest (in my case the Rayeligh).

Is it possible to get such a solution?

I have been able to calculate the first and iterations of the convolution. For the $$N = 1$$ case of the above convolution integral. $$\int_{-\infty}^{+\infty} R(X) R(x - X) \ dX = \\ \frac{1}{4 \sigma^{3}}\exp\left(\frac{-x^{2}}{2\sigma^{2}}\right) \left( 2 x \sigma + \exp\left(\frac{x^{2}}{4\sigma^{2}}\right) \sqrt{\pi} \left( x^{2} - 2 \sigma^{2} \right) \rm{erf}(x/2 \sigma) \right)$$

Where $$R(x)$$ is the Rayleigh distribution. Plotting the above gives a reasonable result if I compare it to my data:

• If you are trying to show that $\sqrt{N}(\bar{X}-\mu)$ has a normal limiting distribution, then using moment generating functions might be an simpler approach.
– JimB
Apr 22, 2020 at 5:07
• .. or, more generally, characteristic functions (/Fourier transforms). Apr 22, 2020 at 5:53
• @JimB are you Mathematica S.E. JimB? Care to provide an answer? I've not come across this technique. Apr 22, 2020 at 14:11
• I'm not sure what you mean by a "real solution." There are a great many useful ways to specify distributions, including (but not limited to) the CDF, PDF, MGF, CF, and CGF. For analyzing them, various approximations and asymptotic expressions can be quite useful. Often it's fruitless to insist on obtaining one version (such as the PDF) but some other version can be extremely simple.
– whuber
Apr 22, 2020 at 16:29
• You have a small typo in the last equation: $\exp{(-x^2/4\sigma^2)}$ should be $\exp{(x^2/4\sigma^2)}$. And that gives you the density for the sum rather than the mean of two i.i.d. Rayleigh random variables. In short, I think you should follow the advice of the two experts. If your objective is to obtain the limiting distribution, then using characteristic functions is certainly the better way to go. If you think you need something close to "exact" distribution for say $N=43$, then comparing the distribution of a huge random sample against the limiting distribution might be adequate.
– JimB
Apr 22, 2020 at 17:25

Here is a quick outline determining the limiting distribution of the sample mean of $$X_1,X_2,\ldots,X_n$$ as $$n\rightarrow \infty$$ using moment generating functions.

All of the $$X_i$$ are independent and have the same distribution: $$X_i \sim \text{Rayleigh}(\sigma)$$. We have $$\mu=\text{E}(X_i)=\sigma \sqrt{\pi/2}$$. Let

$$S_n=\sum_{i=1}^n X_i$$

We look at the limiting distribution of

$$Z=\sqrt{n}(S_n/n-\mu)$$

using moment generating functions. The moment generating function of a Rayleigh distribution is

$$M(t)=\sqrt{\frac{\pi }{2}} \sigma t \left(\text{erf}\left(\frac{\sigma t}{\sqrt{2}}\right)+1\right) \exp \left(\frac{\sigma ^2 t^2}{2}\right)+1$$

The moment generating function of $$S_n$$ is simply $$M(t)^n$$. We also know that the moment generating function of $$a S_n+b$$ (where $$a$$ and $$b$$ are known constants) is $$e^{b t} M(a t)^n$$. So the moment generating function of $$Z=\sqrt{n}(S_n/n-\mu)$$ is

$$e^{\sqrt{\frac{\pi }{2}} \left(-\sqrt{n}\right) \sigma t} \left(\frac{\sqrt{\frac{\pi }{2}} \sigma t e^{\frac{\sigma ^2 t^2}{2 n}} \left(\text{erf}\left(\frac{\sigma t}{\sqrt{2} \sqrt{n}}\right)+1\right)}{\sqrt{n}}+1\right)^n$$

Taking the limit as $$n->\infty$$ we have

$$e^{\frac{1}{4} (4-\pi) \sigma ^2 t^2}$$

for the moment generating function. As the moment generating function for a normal distrbution with mean $$m$$ and variance $$v$$ is

$$e^{m t+\frac{t^2 v}{2}}$$

we see that our limiting moment generating function represents a normal distribution with mean $$0$$ and variance $$\frac{1}{2} (4-\pi ) \sigma ^2$$. What this means is that we can use a normal distribution with mean $$\sigma \sqrt{\pi/2}$$ and variance $$\frac{1}{2} (4-\pi ) \sigma ^2/n$$ to approximate the distribution of $$\bar{X}=S_n/n$$.

As a partial check consider $$n=10$$. We can generate a bunch of samples and plot the resulting histogram and approximate density (using Mathematica)

n = 10;
\[Sigma] = 1;
z = RandomVariate[RayleighDistribution[1], {1000000, n}];
z = Mean[#] & /@ z;
Show[Histogram[z, 100, "PDF", Frame -> True,
FrameLabel -> (Style[#, Bold, 18] &) /@ {"\!$$\*OverscriptBox[\(X$$, $$_$$]\)", "Probability density"}],
Plot[PDF[NormalDistribution[Sqrt[\[Pi]/2] \[Sigma], Sqrt[(4 - \[Pi]) \[Sigma]^2/(2 n)]], x],
{x, 0.4, 3}, PlotRange -> All]]


• Because the OP is aware of the CLT, it's unnecessary to reproduce the steps in a proof: it suffices to compute the mean and variance of the sum of $n$ iid Rayleigh variables and apply the CLT. But I doubt that's what the OP was looking for in an answer.
– whuber
Apr 22, 2020 at 22:47
• @whuber I understand. But the OP (in a comment) mentioned "I've not come across this technique." So I gave a few details.
– JimB
Apr 22, 2020 at 23:55
• I'm going to accept @JimB's answer. This has provided me a lot more information to play with. Thanks! Apr 23, 2020 at 23:12
• If the need is to find pdf values for a particular value of $N$ (given that there doesn't appear to be a nice closed-form solution for $N>2$), then are a few ways to do that. Here are three: (1) Take a large random sample and obtain a nonparmetric density estimate (SmoothHistogram and SmoothKernelDistribution). (2) Use another distribution that matches say the first 3 moments and use that as an approximation. The SkewNormal distribution might be considered (en.wikipedia.org/wiki/Skew_normal_distribution). (3) Use numerical integration (but only useful for small $N$).
– JimB
Apr 24, 2020 at 0:24