# Sampling the max among $N$ samples from the Rayleigh distribution

I've read on the internet that the pdf of the sample max, $$X$$, from among $$N$$ i.i.d. samples from a distribution with pdf $$f(x)$$ and cdf $$F(x)$$ is given by

$$p(X) = N f(x) F(x)^{N-1}.$$

I'm interested in the distribution of the max among $$N$$ samples from a Rayleigh distribution, and I've found numerically that the equation above seems to fit very well, as shown in the figure below.

My question is: Is there an easy and efficient way to sample directly from this distribution? I have been able to obtain samples by rejection sampling, but my implementation is not super fast, so I thought I would ask if there are any other methods I should be aware of. For example, is it known that the distribution of the sample max from a Rayleigh distribution forms some other, well-known distribution?

In case it is relevant: In my application I need to draw very many such samples, hence the need for efficiency, and $$N$$ will vary from 0 to maybe a few hundred.

To quote from Wikipedia

The probability density function of the Rayleigh distribution is$$f(x;\sigma) = \frac{x}{\sigma^2} e^{-x^2/(2\sigma^2)}, \quad x \geq 0,$$where $$\sigma$$ is the scale parameter of the distribution. The cumulative distribution function is$$F(x;\sigma) = 1 - e^{-x^2/(2\sigma^2)}$$for $$x \in (0,\infty).$$

Therefore the density of $$X_{(N)}$$, the largest order statistic, is $$f_{(N)}(x)=N\frac{x}{\sigma^2} e^{-x^2/(2\sigma^2)}[1 - e^{-x^2/(2\sigma^2)}]^{N-1}\quad x\ge 0\tag{1}$$ Simulating from $$f_{(N)}$$ can proceed by

1. acceptance-rejection for (1) if a tight upper bound can be found (see 3.)
2. direct simulation of a sample of $$N$$ Rayleigh variates and derivation of the largest value
3. an expansion of $$[1 - e^{-x^2/(2\sigma^2)}]^{N-1}$$ found in (1) into$$\sum_{i=0}^{N-1} (-1)^i{N-1 \choose i}e^{-ix^2/(2\sigma^2)}$$ and the exploitation of the resulting density as a signed mixture, whose positive elements can be used as proposal in an acceptance-rejection algorithm
4. taking advantage of the representation of a Rayleigh random variate as$$^*$$$$X_i=\sigma\sqrt{-2\log U_i}\qquad U_i\sim\mathcal U(0,1)$$since it implies that the largest value$$X_{(N)}=\sigma\sqrt{-2\log U_{(1)}}\qquad U_{(1)}\sim\mathcal Be(1,N)$$corresponds to the smallest Uniform variate, $$U_{(1)}$$, distributed from the $$\mathcal Be(1,N)$$ Beta distribution.
5. inverting the cdf of $$X_{(N)}$$,$$F_{(N)}(x)=[1 - e^{-x^2/(2\sigma^2)}]^{N}\quad x\ge 0$$into $$X_{(N)}=F_{(N)}^{-1}(U)=\sigma\sqrt{-2\log(1-U^{1/N})}\qquad U\sim\mathcal U(0,1)$$

Options 4 and 5 are equivalent and obviously the most efficient options since their computing time does not increase with $$N$$ and does not involve rejection.

A quick benchmark

  "beta"={
N=339
F=0;for(t in 1:1e3)F=F+rbeta(1,1,N)
},
"inverse"={
N=339;ooN=1/N
F=0;for(t in 1:1e3)F=F+sqrt(-log(1-runif(1)**ooN))
F=F*sqrt(2)
},
"inverse+"={
N=339
F=0;for(t in 1:1e3)F=F+sqrt(-log(1-exp(-rexp(1,N))))
F=F*sqrt(2)
},
replications=1e3}


gives the Beta simulation as fastest (in R):

      test replications elapsed relative user.self sys.self
1     beta         1000   4.208     1.00     4.203    0.003
2  inverse         1000   5.091     1.21     5.092    0.000
3 inverse+         1000   5.304     1.26     5.306    0.000


However, a vectorial version

  "beta"={
N=339
mean(rbeta(1e3,1,N))
},
"inverse"={
N=339;ooN=1/N
mean(sqrt(-log(1-runif(1e3)**ooN)))*sqrt(2)
},
"inverse+"={
N=339
mean(sqrt(-log(1-exp(-rexp(1e3,N)))))*sqrt(2)
}


      test replications elapsed relative user.self sys.self
1     beta         1000   0.142    1.797     0.143        0
2  inverse         1000   0.081    1.025     0.081        0
3 inverse+         1000   0.079    1.000     0.079        0


$$^*$$ This representation follows from the interpretation of the standard Rayleigh distribution as the length of a standard Normal bidimensional vector, $$X=\sqrt{X_1^2+X_2^2}$$, which can be equally written as$$X=\sqrt{-2\cos(2\pi V)^2\log(U)-2\sin(2\pi V)^2\log(U)}$$ as in the Box-Muller algorithm.

• Could you expand a little on point 4? What is $\mathcal{B}e$? Only thing I can think of is the Bernoulli distribution, but then the outcome is either 0 or 1, isn't it?
– Tor
Commented Nov 25, 2023 at 9:21
• That worked like a charm, many thanks! Out of curiosity, where can one learn these tricks? Are there books that explain the relation between different distributions, and how they can be used in practice to draw samples?
– Tor
Commented Nov 25, 2023 at 9:49
• The ultimate reference in simulation remains the 1987 (!) book of Luc Devroye, Non-Uniform Random Variate Generation. Commented Nov 25, 2023 at 14:08
• Quite a few years ago now, Luc made a version of his out of print book available for free on line once the original publisher indicated they had no further interest in it. If you go here: luc.devroye.org/rnbookindex.html you can find links to a pdf of the book and a link to errata for the book. There are newer books (and older books for that matter) but there's much in his book that's worth knowing and it's not always in other references. That said, there's obviously been a few developments in the last 35 years or so, which means it's worth also reading around more widely as well. . . Commented Nov 26, 2023 at 3:27
• . . . This document of his (a book chapter) has some more recent information luc.devroye.org/handbooksimulation1.pdf Commented Nov 26, 2023 at 3:28