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

Parameters:¶

X: scipy.stats.rv_continuous | scipy.stats.rv_histogram¶

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

Y: scipy.stats.rv_continuous | scipy.stats.rv_histogram¶

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

z: ArrayLike | None = None¶

The \(z\) values at which to evaluate the PDF of \(Z\). If None, the \(z\) values are determined based on p.

p: float = 1e-10¶

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 product \(Z = X \cdot 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 = X \cdot Y\) as follows.

The PDF of \(Z\), denoted \(f_Z(z)\), can be derived using the joint distribution of \(X\) and \(Y\). Since \(Z = X \cdot Y\), we express the relationship between \(x\), \(y\), and \(z\) and use a transformation approach.

Let \(Z = X \cdot Y\). The PDF \(f_Z(z)\) is given by

The Jacobian adjustment for this transformation contributes the factor \(\frac{1}{\left| w \right|}\).

Examples

Compute the distribution of the product of two normal random variables.

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

In [2]: Y = scipy.stats.norm(loc=2, scale=1.5)

In [3]: x = np.linspace(-15, 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.multiply_rvs(X, Y).pdf(x), label=r"$X \cdot Y$"); \
   ...: plt.hist(X.rvs(100_000) * Y.rvs(100_000), bins=101, density=True, histtype="step", label=r"$X \cdot Y$ empirical"); \
   ...: plt.legend(); \
   ...: plt.xlabel("Random variable"); \
   ...: plt.ylabel("Probability density"); \
   ...: plt.title("Product of two normal random variables");
   ...: