Improve calculation of random variable arithmetic

src/sdr/_detection/_theory.py

@@ -22,7 +22,7 @@
     verify_not_specified,
     verify_scalar,
 )
-from .._probability import sum_distribution
+from .._probability import add_iid_rvs
 
 
 @export
@@ -244,7 +244,7 @@ def _h0(
             h0 = scipy.stats.norm(0, np.sqrt(sigma2_per))
         elif detector == "linear":
             h0 = scipy.stats.chi(nu, scale=np.sqrt(sigma2_per))
-            h0 = sum_distribution(h0, n_nc)
+            h0 = add_iid_rvs(h0, n_nc)
         elif detector == "square-law":
             h0 = scipy.stats.chi2(nu * n_nc, scale=sigma2_per)
 
@@ -490,7 +490,7 @@ def _h1(
             else:
                 # Folded normal distribution has 1 degree of freedom
                 h1 = scipy.stats.foldnorm(np.sqrt(lambda_), scale=np.sqrt(sigma2_per))
-            h1 = sum_distribution(h1, n_nc)
+            h1 = add_iid_rvs(h1, n_nc)
         elif detector == "square-law":
             h1 = scipy.stats.ncx2(nu * n_nc, lambda_ * n_nc, scale=sigma2_per)
 

src/sdr/_simulation/_channel_model.py

@@ -106,9 +106,9 @@ def dmc(
         x: The input sequence $x$ with $x_i \in \mathcal{X}$.
         P: The $m \times n$ transition probability matrix $P$, where $P_{i,j} = \Pr(Y = y_j | X = x_i)$.
         X: The input alphabet $\mathcal{X}$ of size $m$. If `None`, it is assumed that
-            $\mathcal{X} = \{0, 1, \ldots, m-1\}$.
+            $\mathcal{X} = \{0, 1, \dots, m-1\}$.
         Y: The output alphabet $\mathcal{Y}$ of size $n$. If `None`, it is assumed that
-            $\mathcal{Y} = \{0, 1, \ldots, n-1\}$.
+            $\mathcal{Y} = \{0, 1, \dots, n-1\}$.
         seed: The seed for the random number generator. This is passed to :func:`numpy.random.default_rng`.
 
     Returns:
@@ -546,8 +546,8 @@ def __init__(
 
         Arguments:
             P: The $m \times n$ transition probability matrix $P$, where $P = \Pr(Y = y_j | X = x_i)$.
-            X: The input alphabet $\mathcal{X}$ of size $m$. If `None`, it is assumed that $\mathcal{X} = \{0, 1, \ldots, m-1\}$.
-            Y: The output alphabet $\mathcal{Y}$ of size $n$. If `None`, it is assumed that $\mathcal{Y} = \{0, 1, \ldots, n-1\}$.
+            X: The input alphabet $\mathcal{X}$ of size $m$. If `None`, it is assumed that $\mathcal{X} = \{0, 1, \dots, m-1\}$.
+            Y: The output alphabet $\mathcal{Y}$ of size $n$. If `None`, it is assumed that $\mathcal{Y} = \{0, 1, \dots, n-1\}$.
             seed: The seed for the random number generator. This is passed to :func:`numpy.random.default_rng`.
         """
         P = verify_arraylike(P, float=True, ndim=2, inclusive_min=0, inclusive_max=1)

tests/probability/test_add_iid_rvs.py

@@ -0,0 +1,92 @@
+import numpy as np
+import scipy.stats
+
+import sdr
+
+
+def test_normal():
+    X = scipy.stats.norm(loc=-1, scale=0.5)
+    _verify(X, 2)
+    _verify(X, 3)
+    _verify(X, 4)
+
+
+def test_rayleigh():
+    X = scipy.stats.rayleigh(scale=1)
+    _verify(X, 2)
+    _verify(X, 3)
+    _verify(X, 4)
+
+
+def test_rician():
+    X = scipy.stats.rice(2)
+    _verify(X, 2)
+    _verify(X, 3)
+    _verify(X, 4)
+
+
+def test_exponential():
+    X = scipy.stats.expon(scale=1)
+    _verify(X, 2)
+    _verify(X, 3)
+    _verify(X, 4)
+
+
+def test_gamma():
+    X = scipy.stats.gamma(a=2, scale=1)
+    _verify(X, 2)
+    _verify(X, 3)
+    _verify(X, 4)
+
+
+# def test_poisson():
+#     X = scipy.stats.poisson(mu=2)
+#     _verify(X, 2)
+#     _verify(X, 3)
+#     _verify(X, 4)
+
+
+# def test_bernoulli():
+#     X = scipy.stats.bernoulli(p=0.5)
+#     _verify(X, 2)
+#     _verify(X, 3)
+#     _verify(X, 4)
+
+
+# def test_geomertic():
+#     X = scipy.stats.geom(p=0.5)
+#     _verify(X, 2)
+#     _verify(X, 3)
+#     _verify(X, 4)
+
+
+def test_chi2():
+    X = scipy.stats.chi2(df=2)
+    _verify(X, 2)
+    _verify(X, 3)
+    _verify(X, 4)
+
+
+def _verify(X, n):
+    # Numerically compute the distribution
+    Y = sdr.add_iid_rvs(X, n)
+
+    # Empirically compute the distribution
+    y = X.rvs((100_000, n)).sum(axis=1)
+    hist, bins = np.histogram(y, bins=101, density=True)
+    x = bins[1:] - np.diff(bins) / 2
+
+    if False:
+        import matplotlib.pyplot as plt
+
+        plt.figure()
+        plt.plot(x, X.pdf(x), label="X")
+        plt.plot(x, Y.pdf(x), label="X + ... + X")
+        plt.hist(y, bins=101, cumulative=False, density=True, histtype="step", label="X + .. + X empirical")
+        plt.legend()
+        plt.xlabel("Random variable")
+        plt.ylabel("Probability density")
+        plt.title("Sum of distribution")
+        plt.show()
+
+    assert np.allclose(Y.pdf(x), hist, atol=np.max(hist) * 0.1)

tests/probability/test_sum_distributions.py → tests/probability/test_add_rvs.py

@@ -36,7 +36,7 @@ def test_rayleigh_rician():
 
 def _verify(X, Y):
     # Numerically compute the distribution
-    Z = sdr.sum_distributions(X, Y)
+    Z = sdr.add_rvs(X, Y)
 
     # Empirically compute the distribution
     z = X.rvs(100_000) + Y.rvs(100_000)

tests/probability/test_sum_distribution.py → tests/probability/test_max_iid_rvs.py

@@ -7,30 +7,30 @@
 def test_normal():
     X = scipy.stats.norm(loc=-1, scale=0.5)
     _verify(X, 2)
-    _verify(X, 3)
-    _verify(X, 4)
+    _verify(X, 10)
+    _verify(X, 100)
 
 
 def test_rayleigh():
     X = scipy.stats.rayleigh(scale=1)
     _verify(X, 2)
-    _verify(X, 3)
-    _verify(X, 4)
+    _verify(X, 10)
+    _verify(X, 100)
 
 
 def test_rician():
     X = scipy.stats.rice(2)
     _verify(X, 2)
-    _verify(X, 3)
-    _verify(X, 4)
+    _verify(X, 10)
+    _verify(X, 100)
 
 
 def _verify(X, n):
     # Numerically compute the distribution
-    Y = sdr.sum_distribution(X, n)
+    Y = sdr.max_iid_rvs(X, n)
 
     # Empirically compute the distribution
-    y = X.rvs((100_000, n)).sum(axis=1)
+    y = X.rvs((100_000, n)).max(axis=1)
     hist, bins = np.histogram(y, bins=101, density=True)
     x = bins[1:] - np.diff(bins) / 2
 
@@ -39,8 +39,10 @@ def _verify(X, n):
 
         plt.figure()
         plt.plot(x, X.pdf(x), label="X")
-        plt.plot(x, Y.pdf(x), label="X + ... + X")
-        plt.hist(y, bins=101, cumulative=False, density=True, histtype="step", label="X + .. + X empirical")
+        plt.plot(x, Y.pdf(x), label=r"$\max(X_1, \dots, X_n)$")
+        plt.hist(
+            y, bins=101, cumulative=False, density=True, histtype="step", label=r"$\max(X_1, \dots, X_n)$ empirical"
+        )
         plt.legend()
         plt.xlabel("Random variable")
         plt.ylabel("Probability density")

tests/probability/test_min_iid_rvs.py

@@ -0,0 +1,52 @@
+import numpy as np
+import scipy.stats
+
+import sdr
+
+
+def test_normal():
+    X = scipy.stats.norm(loc=-1, scale=0.5)
+    _verify(X, 2)
+    _verify(X, 10)
+    _verify(X, 100)
+
+
+def test_rayleigh():
+    X = scipy.stats.rayleigh(scale=1)
+    _verify(X, 2)
+    _verify(X, 10)
+    _verify(X, 100)
+
+
+def test_rician():
+    X = scipy.stats.rice(2)
+    _verify(X, 2)
+    _verify(X, 10)
+    _verify(X, 100)
+
+
+def _verify(X, n):
+    # Numerically compute the distribution
+    Y = sdr.min_iid_rvs(X, n)
+
+    # Empirically compute the distribution
+    y = X.rvs((100_000, n)).min(axis=1)
+    hist, bins = np.histogram(y, bins=101, density=True)
+    x = bins[1:] - np.diff(bins) / 2
+
+    if False:
+        import matplotlib.pyplot as plt
+
+        plt.figure()
+        plt.plot(x, X.pdf(x), label="X")
+        plt.plot(x, Y.pdf(x), label=r"$\max(X_1, \dots, X_n)$")
+        plt.hist(
+            y, bins=101, cumulative=False, density=True, histtype="step", label=r"$\max(X_1, \dots, X_n)$ empirical"
+        )
+        plt.legend()
+        plt.xlabel("Random variable")
+        plt.ylabel("Probability density")
+        plt.title("Sum of distribution")
+        plt.show()
+
+    assert np.allclose(Y.pdf(x), hist, atol=np.max(hist) * 0.1)

tests/probability/test_multiply_distribution.py → tests/probability/test_multiply_rvs.py

@@ -41,7 +41,7 @@ def _verify(X, Y):
     x = bins[1:] - np.diff(bins) / 2
 
     # Numerically compute the distribution, only do so over the histogram bins (for speed)
-    Z = sdr.multiply_distributions(X, Y, x)
+    Z = sdr.multiply_rvs(X, Y, x)
 
     if False:
         import matplotlib.pyplot as plt

tests/probability/test_subtract_rvs.py

@@ -0,0 +1,60 @@
+import numpy as np
+import scipy.stats
+
+import sdr
+
+
+def test_normal_normal():
+    X = scipy.stats.norm(loc=-1, scale=0.5)
+    Y = scipy.stats.norm(loc=2, scale=1.5)
+    _verify(X, Y)
+
+
+def test_rayleigh_rayleigh():
+    X = scipy.stats.rayleigh(scale=1)
+    Y = scipy.stats.rayleigh(loc=1, scale=2)
+    _verify(X, Y)
+
+
+def test_rician_rician():
+    X = scipy.stats.rice(2)
+    Y = scipy.stats.rice(3)
+    _verify(X, Y)
+
+
+def test_normal_rayleigh():
+    X = scipy.stats.norm(loc=-1, scale=0.5)
+    Y = scipy.stats.rayleigh(loc=2, scale=1.5)
+    _verify(X, Y)
+
+
+def test_rayleigh_rician():
+    X = scipy.stats.rayleigh(scale=1)
+    Y = scipy.stats.rice(3)
+    _verify(X, Y)
+
+
+def _verify(X, Y):
+    # Numerically compute the distribution
+    Z = sdr.subtract_rvs(X, Y)
+
+    # Empirically compute the distribution
+    z = X.rvs(100_000) - Y.rvs(100_000)
+    hist, bins = np.histogram(z, bins=101, density=True)
+    x = bins[1:] - np.diff(bins) / 2
+
+    if False:
+        import matplotlib.pyplot as plt
+
+        plt.figure()
+        plt.plot(x, X.pdf(x), label="X")
+        plt.plot(x, Y.pdf(x), label="Y")
+        plt.plot(x, Z.pdf(x), label="X - Y")
+        plt.hist(z, bins=101, cumulative=False, density=True, histtype="step", label="X - Y empirical")
+        plt.legend()
+        plt.xlabel("Random variable")
+        plt.ylabel("Probability density")
+        plt.title("Difference of two distributions")
+        plt.show()
+
+    assert np.allclose(Z.pdf(x), hist, atol=np.max(hist) * 0.1)