Distribution of the Sum of 2 Rayleigh Distributed Values

Question

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)
library(extraDistr)

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?

Answer 1

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)
library(extraDistr)

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.

Answer 2

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)