# Distribution of the Sum of 2 Rayleigh Distributed Values

## Problem

I have been trying to determine the distribution of the sum of 2 values sourced from a Rayleigh distribution with the parameter $$\sigma$$. Unfortunately, I am unable to match a computer simulated distribution (see below). I have grown quite frustrated with this now and have decided to reach out.

## Approaches

I have attempted two different approaches.

### Joint PDF Approach (Casella and Berger Second Edition, Section 4.3, pg 184)

Given that $$f(x|\sigma)=\frac{x}{\sigma^2}e^{-\frac{x^2}{2\sigma^2}}$$ is the parent distribution and I choose $$u=a+b$$ and $$v=a-b$$, I get:

$$f(u,v)=\left(\frac{\left(\frac{u+v}{2}\right)}{\sigma^2}e^{-\frac{\left(\frac{u+v}{2}\right)^2}{2\sigma^2}}\right)\left(\frac{\left(\frac{u-v}{2}\right)}{\sigma^2}e^{-\frac{\left(\frac{u-v}{2}\right)^2}{2\sigma^2}}\right)\left|-\frac{1}{2}\right|$$

In this case the $$-\frac{1}{2}$$ on the end of the formula is from the Jacobian. To get the distribution of $$u$$, which is the sum of the two numbers in my case, I need to just marginalize over $$v$$.

\begin{align*} f(u) & = \int_{-\infty}^{\infty}\left(\frac{\left(\frac{u+v}{2}\right)}{\sigma^2}e^{-\frac{\left(\frac{u+v}{2}\right)^2}{2\sigma^2}}\right)\left(\frac{\left(\frac{u-v}{2}\right)}{\sigma^2}e^{-\frac{\left(\frac{u-v}{2}\right)^2}{2\sigma^2}}\right)\left|-\frac{1}{2}\right|dv \\ & = \frac{\sqrt{\pi}}{4\sigma^3}e^{-\frac{u^2}{4\sigma^2}}\left(u^2-2\sigma^2\right) \end{align*}

I should say that it is possible (and quite easy with some tricks) to do this by hand, but I verified my work through Mathematica and it checked out.

### Convolution Approach

Starting from the same starting point of the distribution $$f(x|\sigma)$$, the distribution of the sum of two numbers is the convolution of the parent distribution with itself.

\begin{align*} f(u) & = \int_{-\infty}^{\infty}\left(\frac{\left(u-x\right)}{\sigma^2}e^{-\frac{\left(u-x\right)^2}{2\sigma^2}}\right)\left(\frac{\left(x\right)}{\sigma^2}e^{-\frac{\left(x\right)^2}{2\sigma^2}}\right)dx\\ & = \frac{\sqrt{\pi}}{4\sigma^3}e^{-\frac{u^2}{4\sigma^2}}\left(u^2-2\sigma^2\right) \end{align*}

This integral was much harder and I just decided to do this via Mathematica (I have no shame).

This is great, both methods yielded the same functional form. However, when I try to plot this in R, I get the following.

# R Code and Comparison with Simulation

library(tidyverse)

n = 2
sigma = 2
samps = sapply(seq_len(100000), function(i) sum(rrayleigh(n, sigma=sigma)))

f = function(x) (sqrt(pi)/(4*(sigma**3)))*exp(-1*(x**2)/(4*(sigma**2)))*(x**2-(2*(sigma**2)))

ggplot(data.frame(x=samps)) + theme_bw() + theme(panel.grid=element_blank()) +
stat_density(aes(x=x), color='blue', geom='line') +
stat_function(fun=f)


In the plot, the blue curve is from a quick simulation run and the black curve is the function $$f$$ from above. If I did things correctly, they should match.

As you can see, the function $$f$$ goes negative, which does not make sense for a pdf. This is of course evident in the functional form in the $$\left(u^2-2\sigma^2\right)$$ term which goes negative when $$u=\sqrt{2}\sigma$$.

# Despite Plee for Help

Can somebody here help me understand where I went wrong?

• Your convolution integral is all wrong. The support of the Rayleigh pdf is $[0,\infty)$ and so your convolution should run from $0$ to $u$ and not over the entire real line. – Dilip Sarwate Jun 22 '19 at 3:24
• Copied from this link "Abstract: Sums of Rayleigh random variables occur extensively in wireless communications. A closed-form expression does not exist for the sum distribution and consequently, it is often evaluated numerically or approximated. ... Highly accurate, simple closed-form expressions for the sum distributions and densities are presented. These approximations are valid for a wide range of argument values and number of summands." – BruceET Jun 22 '19 at 9:01

This is too long for a comment. Using Mathematica, we can get the density of the sum of 2 iid Rayleigh distributions as:

$$f(x;\sigma) = \frac{\exp\left(-\frac{x^2}{2\sigma^2}\right)\left(2x\sigma+\exp\left(\frac{x^2}{4\sigma^2}\right)\sqrt{\pi}\left(x^2-2\sigma^2\right)\mathrm{Erf}\left(\frac{x}{2\sigma}\right)\right)}{4\sigma^3}\quad x> 0$$

where $$\mathrm{Erf}$$ denotes the error function:

$$\mathrm{Erf}(x) = \frac{1}{\sqrt{\pi}}\int_{-x}^{x}\exp\left(-t^2\right)dt$$

As @BruceET mentions in the comments, the error function is not available in closed form so the density mentioned above is also not in closed form.

Checking using R:

library(tidyverse)

set.seed(142857)
n = 2
sigma = 2
samps = sapply(seq_len(100000), function(i) sum(rrayleigh(n, sigma=sigma)))

# Error function
erf <- function(x) 2 * pnorm(x * sqrt(2)) - 1
# Density function

# Density function
ray_dens <- Vectorize(function(x, sigma) {
if (x > 0) {
(exp(-x^2/(2*sigma^2))*(2*x*sigma+exp(x^2/(4*sigma^2))*sqrt(pi)*(x^2 - 2*sigma^2)*erf(x/(2*sigma))))/(4*sigma^3)
} else {
0
}
})

ggplot(data.frame(x=samps)) + theme_bw() + theme(panel.grid=element_blank()) +
stat_density(aes(x=x), color='blue', geom='line', size = 1.5) +
stat_function(fun=~ray_dens(.x, sigma = sigma), color = "#F07E00", size = 1.5)


The agreement between simulations (blue line) and the theoretical result (orange line) is excellent.

• This is great. Thanks. But can you tell me where I went wrong with my first integral using the method from Casella and Berger? I did change my limits for the convolution approach to 0 to $\infty$ and got your answer. I used $-\infty$ to $\infty$ in the first integral since I was integrating across $v$ which could span from those extremes. – Justace Clutter Jun 22 '19 at 10:01
• And this doesn't disagree with the claim that there is no closed form because the normal CDF doesn't exist in closed form. But it nice to have this density in a more informative format. (+1) – BruceET Jun 22 '19 at 10:01
• @BruceET That's true and I updated my answer to explicitly mention that fact. Thanks. – COOLSerdash Jun 22 '19 at 10:04

Comment continued: It seems there is no closed form expression for the density. So no wonder you're frustrated. Of course, you can get information for specific cases by simulation. Simulated distribution in R for $$\sigma = 1$$ with kernel density estimator:,

set.seed(2019);  m = 10^6
y1 = sqrt(rchisq(m, 2)); y2 = sqrt(rchisq(m,2))
x = y1 + y2
hist(x, prob=T, br=40, col="skyblue2")
lines(density(x), col="red", lwd=2)

summary(x)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
0.06991 1.83712 2.43511 2.50717 3.09896 7.80869
sd(x)
[1] 0.9261972   # aprx SD(X)
sqrt(4 - pi)
[1] 0.9265028   # exact SD(X)


Note: In case the web link is ephemeral, the reference is: J.Hu, N.C, Breslieu, et al. (2005) Accurate simple closed-form approximations to Rayleigh sum distributions and densities, IEEE Communications Letters Vol. 9, Nr. 2 (Feb)

• Simulation will certainly work but one can do numerical convolution (e.g. via a fast Fourier transform); on a univariate problem it's typically relatively fast and accurate. – Glen_b Jun 22 '19 at 9:41
• With a million iterations simulation gives only about 2 to 3 place accuracy for mean, SD, etc. So your method seems better. Mainly, since OP was struggling with the shape of the density curve, I just wanted to finish my 'bad news' comment by showing the KDE. – BruceET Jun 22 '19 at 9:52
• As far as the no closed form, I guess you mean there is no general closed form for the sum of n random pulls. Right? – Justace Clutter Jun 22 '19 at 9:59
• Also, why did you use the chisel distribution above vice the rayleigh? – Justace Clutter Jun 22 '19 at 10:07
• You can write your own functions for the Rayleigh distribution, perhaps as a first step with something like dray <- function(x, sigma=1){x/sigma^2 * exp(-x^2/(2*sigma^2))}; pray <- function(x, sigma=1){1 - exp(-x^2/(2*sigma^2))}; qray <- function(p, sigma=1){sqrt(-log(1-p)*(2*sigma^2))}; rray <- function(n, sigma=1){qray(runif(n), sigma)} – Henry Jun 22 '19 at 10:42