sdr.multiply_distributions(X: scipy.stats.rv_continuous | scipy.stats.rv_histogram, Y: scipy.stats.rv_continuous | scipy.stats.rv_histogram, x: ArrayLike | None = None, p: float = 1e-10) → scipy.stats.rv_histogram

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 first random variable \(X\).

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

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

x: ArrayLike | None = None

The x values at which to evaluate the PDF of the product. If None, the x 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

The PDF of the product of two independent random variables is calculated as follows.

\[ f_{X \cdot Y}(t) = \int_{0}^{\infty} f_X(k) f_Y(t/k) \frac{1}{k} dk - \int_{-\infty}^{0} f_X(k) f_Y(t/k) \frac{1}{k} dk \]

Examples

Compute the distribution of the product of two normal distributions.

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_distributions(X, Y).pdf(x), label="X * Y"); \
   ...: plt.hist(X.rvs(100_000) * Y.rvs(100_000), bins=101, density=True, histtype="step", label="X * Y empirical"); \
   ...: plt.legend(); \
   ...: plt.xlabel("Random variable"); \
   ...: plt.ylabel("Probability density"); \
   ...: plt.title("Product of two Normal distributions");
   ...:
../../_images/sdr_multiply_distributions_1.png