sdr.max_rvs(X: rv_continuous | rv_histogram, Y: rv_continuous | rv_histogram, p: float = 1e-16) → rv_histogram

Numerically calculates the distribution of the maximum of two independent random variables \(X\) and \(Y\).

Parameters:
X: rv_continuous | rv_histogram

The distribution of the random variable \(X\).

Y: rv_continuous | rv_histogram

The distribution of the random variable \(Y\).

p: float = 1e-16

The probability of exceeding the x axis, on either side, for each distribution. This is used to determine the bounds on the x axis for the numerical convolution. Smaller values of \(p\) will result in more accurate analysis, but will require more computation.

Returns:

The distribution of the maximum \(Z = \max(X, Y)\).

Notes

Given two independent random variables \(X\) and \(Y\) with PDFs \(f_X(x)\) and \(f_Y(y)\), and CDFs \(F_X(x)\) and \(F_Y(y)\), we compute the PDF of \(Z = \max(X, Y)\) as follows.

The CDF of \(Z\), denoted \(F_Z(z)\), is \(F_Z(z) = P(Z \leq z)\). Since \(Z = \max(X, Y)\), the event \(Z \leq z\) occurs if and only if both \(X \leq z\) and \(Y \leq z\). Using independence,

\[F_Z(z) = P(X \leq z) \cdot P(Y \leq z) = F_X(z) \cdot F_Y(z) .\]

The PDF of \(Z\), denoted \(f_Z(z)\), is the derivative of \(F_Z(z)\). Therefore, \(f_Z(z) = \frac{d}{dz} F_Z(z)\). Substituting \(F_Z(z) = F_X(z) \cdot F_Y(z)\) yields

\[f_Z(z) = \frac{d}{dz} \big(F_X(z) \cdot F_Y(z)\big)\]
\[f_Z(z) = f_X(z) \cdot F_Y(z) + f_Y(z) \cdot F_X(z) .\]

Therefore, the PDF of \(Z = \max(X, Y)\) is

\[f_Z(z) = f_X(z) \cdot F_Y(z) + f_Y(z) \cdot F_X(z)\]

where \(F_X(z)\) and \(F_Y(z)\) are the CDFs of the original random variables \(X\) and \(Y\), and \(f_X(z)\) and \(f_Y(z)\) are the PDFs of \(X\) and \(Y\).

Examples

Compute the distribution of the maximum of normal and Rayleigh random variables.

In [1]: X = scipy.stats.norm(loc=3, scale=0.5)

In [2]: Y = scipy.stats.rayleigh(loc=1, scale=2)

In [3]: x = np.linspace(0, 10, 1_001)

In [4]: plt.figure(); \
   ...: plt.plot(x, X.pdf(x), label="$X$"); \
   ...: plt.plot(x, Y.pdf(x), label="$Y$"); \
   ...: plt.plot(x, sdr.max_rvs(X, Y).pdf(x), label=r"$\max(X, Y)$"); \
   ...: plt.hist(np.maximum(X.rvs(100_000), Y.rvs(100_000)), bins=101, density=True, histtype="step", label=r"$\max(X, Y)$ empirical"); \
   ...: plt.legend(); \
   ...: plt.xlabel("Random variable"); \
   ...: plt.ylabel("Probability density"); \
   ...: plt.title("Maximum of normal and Rayleigh random variables");
   ...:
../../_images/sdr_max_rvs_1.png

Compute the distribution of the maximum of Rayleigh and Rician random variables.

In [5]: X = scipy.stats.rayleigh(scale=1)

In [6]: Y = scipy.stats.rice(2)

In [7]: x = np.linspace(0, 6, 1_001)

In [8]: plt.figure(); \
   ...: plt.plot(x, X.pdf(x), label="$X$"); \
   ...: plt.plot(x, Y.pdf(x), label="$Y$"); \
   ...: plt.plot(x, sdr.max_rvs(X, Y).pdf(x), label=r"$\max(X, Y)$"); \
   ...: plt.hist(np.maximum(X.rvs(100_000), Y.rvs(100_000)), bins=101, density=True, histtype="step", label=r"$\max(X, Y)$ empirical"); \
   ...: plt.legend(); \
   ...: plt.xlabel("Random variable"); \
   ...: plt.ylabel("Probability density"); \
   ...: plt.title("Maximum of Rayleigh and Rician random variables");
   ...:
../../_images/sdr_max_rvs_2.png

Compute the distribution of the maximum of Rician and Chi-squared random variables.

In [9]: X = scipy.stats.rice(3)

In [10]: Y = scipy.stats.chi2(3)

In [11]: x = np.linspace(0, 20, 1_001)

In [12]: plt.figure(); \
   ....: plt.plot(x, X.pdf(x), label="$X$"); \
   ....: plt.plot(x, Y.pdf(x), label="$Y$"); \
   ....: plt.plot(x, sdr.max_rvs(X, Y).pdf(x), label=r"$\max(X, Y)$"); \
   ....: plt.hist(np.maximum(X.rvs(100_000), Y.rvs(100_000)), bins=101, density=True, histtype="step", label=r"$\max(X, Y)$ empirical"); \
   ....: plt.legend(); \
   ....: plt.xlim(0, 15); \
   ....: plt.xlabel("Random variable"); \
   ....: plt.ylabel("Probability density"); \
   ....: plt.title("Maximum of Rician and Chi-squared random variables");
   ....:
../../_images/sdr_max_rvs_3.png