4
$\begingroup$

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)

enter image description here

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?

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Dilip Sarwate Jun 22 '19 at 3:24
  • 1
    $\begingroup$ 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." $\endgroup$ – BruceET Jun 22 '19 at 9:01
5
$\begingroup$

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)

Rayleighdensity

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Justace Clutter Jun 22 '19 at 10:01
  • 2
    $\begingroup$ 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) $\endgroup$ – BruceET Jun 22 '19 at 10:01
  • $\begingroup$ @BruceET That's true and I updated my answer to explicitly mention that fact. Thanks. $\endgroup$ – COOLSerdash Jun 22 '19 at 10:04
2
$\begingroup$

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)

enter image description here

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)

$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – Glen_b -Reinstate Monica Jun 22 '19 at 9:41
  • 1
    $\begingroup$ 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. $\endgroup$ – BruceET Jun 22 '19 at 9:52
  • $\begingroup$ 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? $\endgroup$ – Justace Clutter Jun 22 '19 at 9:59
  • $\begingroup$ Also, why did you use the chisel distribution above vice the rayleigh? $\endgroup$ – Justace Clutter Jun 22 '19 at 10:07
  • 1
    $\begingroup$ 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)} $\endgroup$ – Henry Jun 22 '19 at 10:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.