From 19ec6f74bbccd4289b72b324a17123ace872898c Mon Sep 17 00:00:00 2001 From: mhostetter Date: Tue, 17 Dec 2024 19:32:12 -0500 Subject: [PATCH 01/16] Fix `overlap` input condition Fixes #444 --- src/sdr/plot/_spectral_estimation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sdr/plot/_spectral_estimation.py b/src/sdr/plot/_spectral_estimation.py index 9e7edac1f..6b0008be9 100644 --- a/src/sdr/plot/_spectral_estimation.py +++ b/src/sdr/plot/_spectral_estimation.py @@ -64,7 +64,7 @@ def periodogram( sample_rate, sample_rate_provided = verify_sample_rate(sample_rate) # verify_isinstance(window, str, optional=True) verify_scalar(length, optional=True, int=True, positive=True) - verify_scalar(overlap, optional=True, int=True, positive=True) + verify_scalar(overlap, optional=True, int=True, non_negative=True) verify_scalar(fft, optional=True, int=True, positive=True) verify_literal(detrend, ["constant", "linear", False]) verify_literal(average, ["mean", "median"]) @@ -174,7 +174,7 @@ def spectrogram( sample_rate, sample_rate_provided = verify_sample_rate(sample_rate) # verify_isinstance(window, str, optional=True) verify_scalar(length, optional=True, int=True, positive=True) - verify_scalar(overlap, optional=True, int=True, positive=True) + verify_scalar(overlap, optional=True, int=True, non_negative=True) verify_scalar(fft, optional=True, int=True, positive=True) verify_literal(detrend, ["constant", "linear", False]) verify_isinstance(ax, plt.Axes, optional=True) From b96770bafdcec7b73dec4cd4f40cecdf43d39ab3 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Wed, 18 Dec 2024 06:07:42 -0500 Subject: [PATCH 02/16] Rearrange API --- docs/api/misc.rst | 10 ++ src/sdr/_probability.py | 362 ++++++++++++++++++++-------------------- 2 files changed, 191 insertions(+), 181 deletions(-) diff --git a/docs/api/misc.rst b/docs/api/misc.rst index dd17e1a57..74a010a4b 100644 --- a/docs/api/misc.rst +++ b/docs/api/misc.rst @@ -6,6 +6,16 @@ Probability .. python-apigen-group:: probability +Independent random variables +---------------------------- + +.. python-apigen-group:: independent-rvs + +IID random variables +-------------------- + +.. python-apigen-group:: iid-rvs + Data manipulation ----------------- diff --git a/src/sdr/_probability.py b/src/sdr/_probability.py index 5e4adbf65..89f8daa06 100644 --- a/src/sdr/_probability.py +++ b/src/sdr/_probability.py @@ -77,182 +77,6 @@ def Qinv(p: npt.ArrayLike) -> npt.NDArray[np.float64]: return convert_output(x) -@export -def add_iid_rvs( - X: scipy.stats.rv_continuous | scipy.stats.rv_histogram, - n_vars: int, - p: float = 1e-16, -) -> scipy.stats.rv_histogram: - r""" - Numerically calculates the distribution of the sum of $n$ i.i.d. random variables $X_i$. - - Arguments: - X: The distribution of the i.i.d. random variables $X_i$. - n_vars: The number $n$ of random variables. - p: 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 sum $Z = X_1 + X_2 + \dots + X_n$. - - Notes: - Given a random variable $X$ with PDF $f_X(x)$, we compute the PDF of - $Z = X_1 + X_2 + \dots + X_n$, where $X_1, X_2, \dots, X_n$ are independent and identically distributed - (i.i.d.), as follows. - - The PDF of $Z$, denoted $f_Z(z)$, can be obtained by using the convolution formula for independent - random variables. Specifically, the PDF of the sum of $n$ i.i.d. random variables is given by the $n$-fold - convolution of the PDF of $X$ with itself. - - For $n = 2$, $Z = X_1 + X_2$. The PDF of $Z$ is - - $$f_Z(z) = \int_{-\infty}^\infty f_X(x) f_X(z - x) \, dx$$ - - For $n > 2$, the PDF of $Z = X_1 + X_2 + \dots + X_n$ is computed recursively - - $$f_Z(z) = \int_{-\infty}^\infty f_X(x) f_{Z_{n-1}}(z - x) \, dx$$ - - where $f_{Z_{n-1}}(z)$ is the PDF of the sum of $n-1$ random variables. - - For large $n$, the Central Limit Theorem may be used as an approximation. If $X_i$ have mean $\mu$ and - variance $\sigma^2$, then $Z$ approximately follows $Z \sim \mathcal{N}(n\mu, n\sigma^2)$ - for sufficiently large $n$. - - Examples: - Compute the distribution of the sum of two normal random variables. - - .. ipython:: python - - X = scipy.stats.norm(loc=-1, scale=0.5) - n_vars = 2 - x = np.linspace(-6, 2, 1_001) - - @savefig sdr_add_iid_rvs_1.png - plt.figure(); \ - plt.plot(x, X.pdf(x), label="X"); \ - plt.plot(x, sdr.add_iid_rvs(X, n_vars).pdf(x), label="$X_1 + X_2$"); \ - plt.hist(X.rvs((100_000, n_vars)).sum(axis=1), bins=101, density=True, histtype="step", label="$X_1 + X_2$ empirical"); \ - plt.legend(); \ - plt.xlabel("Random variable"); \ - plt.ylabel("Probability density"); \ - plt.title("Sum of two normal random variables"); - - Compute the distribution of the sum of three Rayleigh random variables. - - .. ipython:: python - - X = scipy.stats.rayleigh(scale=1) - n_vars = 3 - x = np.linspace(0, 10, 1_001) - - @savefig sdr_add_iid_rvs_2.png - plt.figure(); \ - plt.plot(x, X.pdf(x), label="X"); \ - plt.plot(x, sdr.add_iid_rvs(X, n_vars).pdf(x), label="$X_1 + X_2 + X_3$"); \ - plt.hist(X.rvs((100_000, n_vars)).sum(axis=1), bins=101, density=True, histtype="step", label="$X_1 + X_2 + X_3$ empirical"); \ - plt.legend(); \ - plt.xlabel("Random variable"); \ - plt.ylabel("Probability density"); \ - plt.title("Sum of three Rayleigh random variables"); - - Compute the distribution of the sum of four Rician random variables. - - .. ipython:: python - - X = scipy.stats.rice(2) - n_vars = 4 - x = np.linspace(0, 18, 1_001) - - @savefig sdr_add_iid_rvs_3.png - plt.figure(); \ - plt.plot(x, X.pdf(x), label="X"); \ - plt.plot(x, sdr.add_iid_rvs(X, n_vars).pdf(x), label="$X_1 + X_2 + X_3 + X_4$"); \ - plt.hist(X.rvs((100_000, n_vars)).sum(axis=1), bins=101, density=True, histtype="step", label="$X_1 + X_2 + X_3 + X_4$ empirical"); \ - plt.legend(); \ - plt.xlabel("Random variable"); \ - plt.ylabel("Probability density"); \ - plt.title("Sum of four Rician random variables"); - - Group: - probability - """ - verify_scalar(n_vars, int=True, positive=True) - verify_scalar(p, float=True, exclusive_min=0, inclusive_max=0.1) - - if n_vars == 1: - return X - - # Compute the exact distribution, if possible - if isinstance(X.dist, scipy.stats.rv_continuous): - shape, loc, scale = X.dist._parse_args(*X.args, **X.kwds) - - if isinstance(X.dist, type(scipy.stats.norm)): - # Sum of normals is normal - mu = loc # Mean - sigma = scale # Standard deviation - mu_sum = n_vars * mu - sigma_sum = np.sqrt(n_vars) * sigma - return scipy.stats.norm(loc=mu_sum, scale=sigma_sum) - if isinstance(X.dist, type(scipy.stats.expon)): - # Sum of exponentials is gamma - lambda_inv = 1 / X.mean() # Rate parameter lambda - return scipy.stats.gamma(a=n_vars, scale=1 / lambda_inv) - if isinstance(X.dist, type(scipy.stats.gamma)): - # Sum of gammas with the same scale is gamma - a = shape[0] # Shape parameter a - a_sum = a * n_vars - return scipy.stats.gamma(a=a_sum, scale=scale) - if isinstance(X.dist, type(scipy.stats.chi2)): - # Sum of Chi-Squares is Chi-Square - df = shape[0] - df_sum = n_vars * df - return scipy.stats.chi2(df=df_sum) - # if isinstance(X.dist, scipy.stats.rv_discrete): - # shape, loc, scale = X.dist._parse_args(*X.args, **X.kwds) - - # if isinstance(X.dist, type(scipy.stats.poisson)): - # # Sum of Poissons is Poisson - # mu = shape[0] # Rate parameter mu - # mu_sum = mu * n_vars - # return scipy.stats.poisson(mu=mu_sum) - # if isinstance(X.dist, type(scipy.stats.bernoulli)): - # # Sum of Bernoullis is Binomial - # p = shape[1] # Probability of success - # return scipy.stats.binom(n=n_vars, p=p) - # if isinstance(X.dist, type(scipy.stats.geom)): - # # Sum of Geometrics is Negative Binomial - # p = shape[0] # Probability of success - # return scipy.stats.binom(n=n_vars, p=p) - - # TODO: Add an override that uses the Central Limit Theorem for large values of n - - # Determine the x axis of each distribution such that the probability of exceeding the x axis, on either side, - # is p. - z1_min, z1_max = _x_range(X, p) - z = np.linspace(z1_min, z1_max, 1_001) - dz = np.mean(np.diff(z)) - - # Compute the PDF of the base distribution - f_X = X.pdf(z) - - # The PDF of the sum of n_vars independent random variables is the convolution of the PDF of the base distribution. - # This is efficiently computed in the frequency domain by exponentiating the FFT. The FFT must be zero-padded - # enough that the circular convolutions do not wrap around. - n_fft = scipy.fft.next_fast_len(f_X.size * n_vars) - f_X_fft = np.fft.fft(f_X, n_fft) - f_X_fft = f_X_fft**n_vars - f_Z = np.fft.ifft(f_X_fft).real - f_Z /= f_Z.sum() * dz - z = np.arange(f_Z.size) * dz + z[0] * (n_vars) - - # Adjust the histograms bins to be on either side of each point. So there is one extra point added. - z = np.append(z, z[-1] + dz) - z -= dz / 2 - - return scipy.stats.rv_histogram((f_Z, z)) - - @export def add_rvs( X: scipy.stats.rv_continuous | scipy.stats.rv_histogram, @@ -340,7 +164,7 @@ def add_rvs( plt.title("Sum of Rician and Chi-squared random variables"); Group: - probability + independent-rvs """ verify_scalar(p, float=True, exclusive_min=0, inclusive_max=0.1) @@ -459,7 +283,7 @@ def subtract_rvs( plt.title("Difference of Rician and Chi-squared random variables"); Group: - probability + independent-rvs """ verify_scalar(p, float=True, exclusive_min=0, inclusive_max=0.1) @@ -546,7 +370,7 @@ def multiply_rvs( plt.title("Product of two normal random variables"); Group: - probability + independent-rvs """ verify_scalar(p, float=True, exclusive_min=0, inclusive_max=0.1) @@ -578,6 +402,182 @@ def integrand(y: float, z: float) -> float: return scipy.stats.rv_histogram((f_Z, z)) +@export +def add_iid_rvs( + X: scipy.stats.rv_continuous | scipy.stats.rv_histogram, + n_vars: int, + p: float = 1e-16, +) -> scipy.stats.rv_histogram: + r""" + Numerically calculates the distribution of the sum of $n$ i.i.d. random variables $X_i$. + + Arguments: + X: The distribution of the i.i.d. random variables $X_i$. + n_vars: The number $n$ of random variables. + p: 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 sum $Z = X_1 + X_2 + \dots + X_n$. + + Notes: + Given a random variable $X$ with PDF $f_X(x)$, we compute the PDF of + $Z = X_1 + X_2 + \dots + X_n$, where $X_1, X_2, \dots, X_n$ are independent and identically distributed + (i.i.d.), as follows. + + The PDF of $Z$, denoted $f_Z(z)$, can be obtained by using the convolution formula for independent + random variables. Specifically, the PDF of the sum of $n$ i.i.d. random variables is given by the $n$-fold + convolution of the PDF of $X$ with itself. + + For $n = 2$, $Z = X_1 + X_2$. The PDF of $Z$ is + + $$f_Z(z) = \int_{-\infty}^\infty f_X(x) f_X(z - x) \, dx$$ + + For $n > 2$, the PDF of $Z = X_1 + X_2 + \dots + X_n$ is computed recursively + + $$f_Z(z) = \int_{-\infty}^\infty f_X(x) f_{Z_{n-1}}(z - x) \, dx$$ + + where $f_{Z_{n-1}}(z)$ is the PDF of the sum of $n-1$ random variables. + + For large $n$, the Central Limit Theorem may be used as an approximation. If $X_i$ have mean $\mu$ and + variance $\sigma^2$, then $Z$ approximately follows $Z \sim \mathcal{N}(n\mu, n\sigma^2)$ + for sufficiently large $n$. + + Examples: + Compute the distribution of the sum of two normal random variables. + + .. ipython:: python + + X = scipy.stats.norm(loc=-1, scale=0.5) + n_vars = 2 + x = np.linspace(-6, 2, 1_001) + + @savefig sdr_add_iid_rvs_1.png + plt.figure(); \ + plt.plot(x, X.pdf(x), label="X"); \ + plt.plot(x, sdr.add_iid_rvs(X, n_vars).pdf(x), label="$X_1 + X_2$"); \ + plt.hist(X.rvs((100_000, n_vars)).sum(axis=1), bins=101, density=True, histtype="step", label="$X_1 + X_2$ empirical"); \ + plt.legend(); \ + plt.xlabel("Random variable"); \ + plt.ylabel("Probability density"); \ + plt.title("Sum of two normal random variables"); + + Compute the distribution of the sum of three Rayleigh random variables. + + .. ipython:: python + + X = scipy.stats.rayleigh(scale=1) + n_vars = 3 + x = np.linspace(0, 10, 1_001) + + @savefig sdr_add_iid_rvs_2.png + plt.figure(); \ + plt.plot(x, X.pdf(x), label="X"); \ + plt.plot(x, sdr.add_iid_rvs(X, n_vars).pdf(x), label="$X_1 + X_2 + X_3$"); \ + plt.hist(X.rvs((100_000, n_vars)).sum(axis=1), bins=101, density=True, histtype="step", label="$X_1 + X_2 + X_3$ empirical"); \ + plt.legend(); \ + plt.xlabel("Random variable"); \ + plt.ylabel("Probability density"); \ + plt.title("Sum of three Rayleigh random variables"); + + Compute the distribution of the sum of four Rician random variables. + + .. ipython:: python + + X = scipy.stats.rice(2) + n_vars = 4 + x = np.linspace(0, 18, 1_001) + + @savefig sdr_add_iid_rvs_3.png + plt.figure(); \ + plt.plot(x, X.pdf(x), label="X"); \ + plt.plot(x, sdr.add_iid_rvs(X, n_vars).pdf(x), label="$X_1 + X_2 + X_3 + X_4$"); \ + plt.hist(X.rvs((100_000, n_vars)).sum(axis=1), bins=101, density=True, histtype="step", label="$X_1 + X_2 + X_3 + X_4$ empirical"); \ + plt.legend(); \ + plt.xlabel("Random variable"); \ + plt.ylabel("Probability density"); \ + plt.title("Sum of four Rician random variables"); + + Group: + iid-rvs + """ + verify_scalar(n_vars, int=True, positive=True) + verify_scalar(p, float=True, exclusive_min=0, inclusive_max=0.1) + + if n_vars == 1: + return X + + # Compute the exact distribution, if possible + if isinstance(X.dist, scipy.stats.rv_continuous): + shape, loc, scale = X.dist._parse_args(*X.args, **X.kwds) + + if isinstance(X.dist, type(scipy.stats.norm)): + # Sum of normals is normal + mu = loc # Mean + sigma = scale # Standard deviation + mu_sum = n_vars * mu + sigma_sum = np.sqrt(n_vars) * sigma + return scipy.stats.norm(loc=mu_sum, scale=sigma_sum) + if isinstance(X.dist, type(scipy.stats.expon)): + # Sum of exponentials is gamma + lambda_inv = 1 / X.mean() # Rate parameter lambda + return scipy.stats.gamma(a=n_vars, scale=1 / lambda_inv) + if isinstance(X.dist, type(scipy.stats.gamma)): + # Sum of gammas with the same scale is gamma + a = shape[0] # Shape parameter a + a_sum = a * n_vars + return scipy.stats.gamma(a=a_sum, scale=scale) + if isinstance(X.dist, type(scipy.stats.chi2)): + # Sum of Chi-Squares is Chi-Square + df = shape[0] + df_sum = n_vars * df + return scipy.stats.chi2(df=df_sum) + # if isinstance(X.dist, scipy.stats.rv_discrete): + # shape, loc, scale = X.dist._parse_args(*X.args, **X.kwds) + + # if isinstance(X.dist, type(scipy.stats.poisson)): + # # Sum of Poissons is Poisson + # mu = shape[0] # Rate parameter mu + # mu_sum = mu * n_vars + # return scipy.stats.poisson(mu=mu_sum) + # if isinstance(X.dist, type(scipy.stats.bernoulli)): + # # Sum of Bernoullis is Binomial + # p = shape[1] # Probability of success + # return scipy.stats.binom(n=n_vars, p=p) + # if isinstance(X.dist, type(scipy.stats.geom)): + # # Sum of Geometrics is Negative Binomial + # p = shape[0] # Probability of success + # return scipy.stats.binom(n=n_vars, p=p) + + # TODO: Add an override that uses the Central Limit Theorem for large values of n + + # Determine the x axis of each distribution such that the probability of exceeding the x axis, on either side, + # is p. + z1_min, z1_max = _x_range(X, p) + z = np.linspace(z1_min, z1_max, 1_001) + dz = np.mean(np.diff(z)) + + # Compute the PDF of the base distribution + f_X = X.pdf(z) + + # The PDF of the sum of n_vars independent random variables is the convolution of the PDF of the base distribution. + # This is efficiently computed in the frequency domain by exponentiating the FFT. The FFT must be zero-padded + # enough that the circular convolutions do not wrap around. + n_fft = scipy.fft.next_fast_len(f_X.size * n_vars) + f_X_fft = np.fft.fft(f_X, n_fft) + f_X_fft = f_X_fft**n_vars + f_Z = np.fft.ifft(f_X_fft).real + f_Z /= f_Z.sum() * dz + z = np.arange(f_Z.size) * dz + z[0] * (n_vars) + + # Adjust the histograms bins to be on either side of each point. So there is one extra point added. + z = np.append(z, z[-1] + dz) + z -= dz / 2 + + return scipy.stats.rv_histogram((f_Z, z)) + + @export def min_iid_rvs( X: scipy.stats.rv_continuous | scipy.stats.rv_histogram, @@ -673,7 +673,7 @@ def min_iid_rvs( plt.title("Minimum of 100 Rician random variables"); Group: - probability + iid-rvs """ verify_scalar(n_vars, int=True, positive=True) verify_scalar(p, float=True, exclusive_min=0, inclusive_max=0.1) @@ -795,7 +795,7 @@ def max_iid_rvs( plt.title("Maximum of 100 Rician random variables"); Group: - probability + iid-rvs """ verify_scalar(n_vars, int=True, positive=True) verify_scalar(p, float=True, exclusive_min=0, inclusive_max=0.1) From 16b7298effa47d7c38846ad19b7855e5e78e683b Mon Sep 17 00:00:00 2001 From: mhostetter Date: Wed, 18 Dec 2024 06:29:57 -0500 Subject: [PATCH 03/16] Make labels use math equations --- src/sdr/_probability.py | 46 ++++++++++++++++++++--------------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/src/sdr/_probability.py b/src/sdr/_probability.py index 89f8daa06..7bb8a0c9c 100644 --- a/src/sdr/_probability.py +++ b/src/sdr/_probability.py @@ -116,8 +116,8 @@ def add_rvs( @savefig sdr_add_rvs_1.png plt.figure(); \ - plt.plot(x, X.pdf(x), label="X"); \ - plt.plot(x, Y.pdf(x), label="Y"); \ + plt.plot(x, X.pdf(x), label="$X$"); \ + plt.plot(x, Y.pdf(x), label="$Y$"); \ plt.plot(x, sdr.add_rvs(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(); \ @@ -135,8 +135,8 @@ def add_rvs( @savefig sdr_add_rvs_2.png plt.figure(); \ - plt.plot(x, X.pdf(x), label="X"); \ - plt.plot(x, Y.pdf(x), label="Y"); \ + plt.plot(x, X.pdf(x), label="$X$"); \ + plt.plot(x, Y.pdf(x), label="$Y$"); \ plt.plot(x, sdr.add_rvs(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(); \ @@ -154,8 +154,8 @@ def add_rvs( @savefig sdr_add_rvs_3.png plt.figure(); \ - plt.plot(x, X.pdf(x), label="X"); \ - plt.plot(x, Y.pdf(x), label="Y"); \ + plt.plot(x, X.pdf(x), label="$X$"); \ + plt.plot(x, Y.pdf(x), label="$Y$"); \ plt.plot(x, sdr.add_rvs(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(); \ @@ -234,8 +234,8 @@ def subtract_rvs( @savefig sdr_subtract_rvs_1.png plt.figure(); \ - plt.plot(x, X.pdf(x), label="X"); \ - plt.plot(x, Y.pdf(x), label="Y"); \ + plt.plot(x, X.pdf(x), label="$X$"); \ + plt.plot(x, Y.pdf(x), label="$Y$"); \ plt.plot(x, sdr.subtract_rvs(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(); \ @@ -253,8 +253,8 @@ def subtract_rvs( @savefig sdr_subtract_rvs_2.png plt.figure(); \ - plt.plot(x, X.pdf(x), label="X"); \ - plt.plot(x, Y.pdf(x), label="Y"); \ + plt.plot(x, X.pdf(x), label="$X$"); \ + plt.plot(x, Y.pdf(x), label="$Y$"); \ plt.plot(x, sdr.subtract_rvs(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(); \ @@ -272,8 +272,8 @@ def subtract_rvs( @savefig sdr_subtract_rvs_3.png plt.figure(); \ - plt.plot(x, X.pdf(x), label="X"); \ - plt.plot(x, Y.pdf(x), label="Y"); \ + plt.plot(x, X.pdf(x), label="$X$"); \ + plt.plot(x, Y.pdf(x), label="$Y$"); \ plt.plot(x, sdr.subtract_rvs(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(); \ @@ -360,8 +360,8 @@ def multiply_rvs( @savefig sdr_multiply_rvs_1.png plt.figure(); \ - plt.plot(x, X.pdf(x), label="X"); \ - plt.plot(x, Y.pdf(x), label="Y"); \ + 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(); \ @@ -455,7 +455,7 @@ def add_iid_rvs( @savefig sdr_add_iid_rvs_1.png plt.figure(); \ - plt.plot(x, X.pdf(x), label="X"); \ + plt.plot(x, X.pdf(x), label="$X$"); \ plt.plot(x, sdr.add_iid_rvs(X, n_vars).pdf(x), label="$X_1 + X_2$"); \ plt.hist(X.rvs((100_000, n_vars)).sum(axis=1), bins=101, density=True, histtype="step", label="$X_1 + X_2$ empirical"); \ plt.legend(); \ @@ -473,7 +473,7 @@ def add_iid_rvs( @savefig sdr_add_iid_rvs_2.png plt.figure(); \ - plt.plot(x, X.pdf(x), label="X"); \ + plt.plot(x, X.pdf(x), label="$X$"); \ plt.plot(x, sdr.add_iid_rvs(X, n_vars).pdf(x), label="$X_1 + X_2 + X_3$"); \ plt.hist(X.rvs((100_000, n_vars)).sum(axis=1), bins=101, density=True, histtype="step", label="$X_1 + X_2 + X_3$ empirical"); \ plt.legend(); \ @@ -491,7 +491,7 @@ def add_iid_rvs( @savefig sdr_add_iid_rvs_3.png plt.figure(); \ - plt.plot(x, X.pdf(x), label="X"); \ + plt.plot(x, X.pdf(x), label="$X$"); \ plt.plot(x, sdr.add_iid_rvs(X, n_vars).pdf(x), label="$X_1 + X_2 + X_3 + X_4$"); \ plt.hist(X.rvs((100_000, n_vars)).sum(axis=1), bins=101, density=True, histtype="step", label="$X_1 + X_2 + X_3 + X_4$ empirical"); \ plt.legend(); \ @@ -628,7 +628,7 @@ def min_iid_rvs( @savefig sdr_min_iid_rvs_1.png plt.figure(); \ - plt.plot(x, X.pdf(x), label="X"); \ + plt.plot(x, X.pdf(x), label="$X$"); \ plt.plot(x, sdr.min_iid_rvs(X, n_vars).pdf(x), label=r"$\min(X_1, X_2)$"); \ plt.hist(X.rvs((100_000, n_vars)).min(axis=1), bins=101, density=True, histtype="step", label=r"$\min(X_1, X_2)$ empirical"); \ plt.legend(); \ @@ -646,7 +646,7 @@ def min_iid_rvs( @savefig sdr_min_iid_rvs_2.png plt.figure(); \ - plt.plot(x, X.pdf(x), label="X"); \ + plt.plot(x, X.pdf(x), label="$X$"); \ plt.plot(x, sdr.min_iid_rvs(X, n_vars).pdf(x), label="$\\min(X_1, \dots, X_3)$"); \ plt.hist(X.rvs((100_000, n_vars)).min(axis=1), bins=101, density=True, histtype="step", label="$\\min(X_1, \dots, X_3)$ empirical"); \ plt.legend(); \ @@ -664,7 +664,7 @@ def min_iid_rvs( @savefig sdr_min_iid_rvs_3.png plt.figure(); \ - plt.plot(x, X.pdf(x), label="X"); \ + plt.plot(x, X.pdf(x), label="$X$"); \ plt.plot(x, sdr.min_iid_rvs(X, n_vars).pdf(x), label=r"$\min(X_1, \dots, X_{100})$"); \ plt.hist(X.rvs((100_000, n_vars)).min(axis=1), bins=101, density=True, histtype="step", label=r"$\min(X_1, \dots, X_{100})$ empirical"); \ plt.legend(); \ @@ -750,7 +750,7 @@ def max_iid_rvs( @savefig sdr_max_iid_rvs_1.png plt.figure(); \ - plt.plot(x, X.pdf(x), label="X"); \ + plt.plot(x, X.pdf(x), label="$X$"); \ plt.plot(x, sdr.max_iid_rvs(X, n_vars).pdf(x), label=r"$\max(X_1, X_2)$"); \ plt.hist(X.rvs((100_000, n_vars)).max(axis=1), bins=101, density=True, histtype="step", label=r"$\max(X_1, X_2)$ empirical"); \ plt.legend(); \ @@ -768,7 +768,7 @@ def max_iid_rvs( @savefig sdr_max_iid_rvs_2.png plt.figure(); \ - plt.plot(x, X.pdf(x), label="X"); \ + plt.plot(x, X.pdf(x), label="$X$"); \ plt.plot(x, sdr.max_iid_rvs(X, n_vars).pdf(x), label="$\\max(X_1, \dots, X_3)$"); \ plt.hist(X.rvs((100_000, n_vars)).max(axis=1), bins=101, density=True, histtype="step", label="$\\max(X_1, \dots, X_3)$ empirical"); \ plt.legend(); \ @@ -786,7 +786,7 @@ def max_iid_rvs( @savefig sdr_max_iid_rvs_3.png plt.figure(); \ - plt.plot(x, X.pdf(x), label="X"); \ + plt.plot(x, X.pdf(x), label="$X$"); \ plt.plot(x, sdr.max_iid_rvs(X, n_vars).pdf(x), label=r"$\max(X_1, \dots, X_{100})$"); \ plt.hist(X.rvs((100_000, n_vars)).max(axis=1), bins=101, density=True, histtype="step", label=r"$\max(X_1, \dots, X_{100})$ empirical"); \ plt.legend(); \ From 9374db9b9a547bf2b3fc8d8f945c7085a452298b Mon Sep 17 00:00:00 2001 From: mhostetter Date: Wed, 18 Dec 2024 06:43:07 -0500 Subject: [PATCH 04/16] Add `sdr.max_rvs()` Fixes #443 --- src/sdr/_probability.py | 134 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 134 insertions(+) diff --git a/src/sdr/_probability.py b/src/sdr/_probability.py index 7bb8a0c9c..ee1b45f55 100644 --- a/src/sdr/_probability.py +++ b/src/sdr/_probability.py @@ -402,6 +402,140 @@ def integrand(y: float, z: float) -> float: return scipy.stats.rv_histogram((f_Z, z)) +@export +def max_rvs( + X: scipy.stats.rv_continuous | scipy.stats.rv_histogram, + Y: scipy.stats.rv_continuous | scipy.stats.rv_histogram, + p: float = 1e-16, +) -> scipy.stats.rv_histogram: + r""" + Numerically calculates the distribution of the maximum of two independent random variables $X$ and $Y$. + + Arguments: + X: The distribution of the random variable $X$. + Y: The distribution of the random variable $Y$. + p: 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. + + .. ipython:: python + + X = scipy.stats.norm(loc=3, scale=0.5) + Y = scipy.stats.rayleigh(loc=1, scale=2) + x = np.linspace(0, 10, 1_001) + + @savefig sdr_max_rvs_1.png + 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"); + + Compute the distribution of the maximum of Rayleigh and Rician random variables. + + .. ipython:: python + + X = scipy.stats.rayleigh(scale=1) + Y = scipy.stats.rice(2) + x = np.linspace(0, 6, 1_001) + + @savefig sdr_max_rvs_2.png + 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"); + + Compute the distribution of the maximum of Rician and Chi-squared random variables. + + .. ipython:: python + + X = scipy.stats.rice(3) + Y = scipy.stats.chi2(3) + x = np.linspace(0, 20, 1_001) + + @savefig sdr_max_rvs_2.png + 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"); + + Group: + independent-rvs + """ + verify_scalar(p, float=True, exclusive_min=0, inclusive_max=0.1) + + # Determine the x axis of each distribution such that the probability of exceeding the x axis, on either side, + # is p. + x_min, x_max = _x_range(X, p) + y_min, y_max = _x_range(Y, p) + dx = (x_max - x_min) / 1_000 + dy = (y_max - y_min) / 1_000 + dz = np.min([dx, dy]) # Use the smaller delta x -- must use the same dx for both distributions + z = np.arange(np.min([x_min, y_min]), np.max([x_max, y_max]), dz) + + # Compute the PDF of each distribution + f_X = X.pdf(z) + f_Y = Y.pdf(z) + + # Compute the CDF of each distribution + F_X = X.cdf(z) + F_Y = Y.cdf(z) + + # The PDF of the min of X and Y + f_Z = f_X * F_Y + f_Y * F_X + + # Adjust the histograms bins to be on either side of each point. So there is one extra point added. + z = np.append(z, z[-1] + dz) + z -= dz / 2 + + return scipy.stats.rv_histogram((f_Z, z)) + + @export def add_iid_rvs( X: scipy.stats.rv_continuous | scipy.stats.rv_histogram, From c14a56856d82f7c6ea0c99c3835495dcbeb95978 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Wed, 18 Dec 2024 06:43:19 -0500 Subject: [PATCH 05/16] Add unit tests for #443 --- tests/probability/test_max_rvs.py | 60 +++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) create mode 100644 tests/probability/test_max_rvs.py diff --git a/tests/probability/test_max_rvs.py b/tests/probability/test_max_rvs.py new file mode 100644 index 000000000..fd14e1e4f --- /dev/null +++ b/tests/probability/test_max_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.max_rvs(X, Y) + + # Empirically compute the distribution + z = np.maximum(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("Sum of two distributions") + plt.show() + + assert np.allclose(Z.pdf(x), hist, atol=np.max(hist) * 0.1) From 4d2f4bdc1ba8ebcc0683a4cb58a714645be325a6 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Wed, 18 Dec 2024 06:58:21 -0500 Subject: [PATCH 06/16] Add `sdr.min_rvs()` Fixes #442 --- src/sdr/_probability.py | 139 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 139 insertions(+) diff --git a/src/sdr/_probability.py b/src/sdr/_probability.py index ee1b45f55..308777832 100644 --- a/src/sdr/_probability.py +++ b/src/sdr/_probability.py @@ -402,6 +402,145 @@ def integrand(y: float, z: float) -> float: return scipy.stats.rv_histogram((f_Z, z)) +@export +def min_rvs( + X: scipy.stats.rv_continuous | scipy.stats.rv_histogram, + Y: scipy.stats.rv_continuous | scipy.stats.rv_histogram, + p: float = 1e-16, +) -> scipy.stats.rv_histogram: + r""" + Numerically calculates the distribution of the minimum of two independent random variables $X$ and $Y$. + + Arguments: + X: The distribution of the random variable $X$. + Y: The distribution of the random variable $Y$. + p: 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 minimum $Z = \min(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 = \min(X, Y)$ as follows. + + The CDF of $Z$, denoted $F_Z(z)$, is $F_Z(z) = P(Z \leq z)$. Since $Z = \min(X, Y)$, the event $Z \leq z$ + occurs if either $X \leq z$ or $Y \leq z$. Using the complement and independence, + + $$ + P(Z \leq z) = 1 - P(Z > z) = 1 - P(X > z \text{ and } Y > z) + F_Z(z) = 1 - P(X > z) \cdot P(Y > z) + F_Z(z) = 1 - (1 - F_X(z)) \cdot (1 - 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) = 1 - (1 - F_X(z)) \cdot (1 - F_Y(z))$ yields + + $$ + f_Z(z) = \frac{d}{dz} \big(1 - (1 - F_X(z)) \cdot (1 - F_Y(z))\big) + f_Z(z) = \frac{d}{dz} \big((1 - F_X(z)) \cdot (1 - F_Y(z))\big) + f_Z(z) = f_X(z) \cdot (1 - F_Y(z)) + f_Y(z) \cdot (1 - F_X(z)) . + $$ + + Therefore, the PDF of $Z = \min(X, Y)$ is + + $$f_Z(z) = f_X(z) \cdot (1 - F_Y(z)) + f_Y(z) \cdot (1 - 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 minimum of normal and Rayleigh random variables. + + .. ipython:: python + + X = scipy.stats.norm(loc=3, scale=0.5) + Y = scipy.stats.rayleigh(loc=1, scale=2) + x = np.linspace(0, 10, 1_001) + + @savefig sdr_min_rvs_1.png + plt.figure(); \ + plt.plot(x, X.pdf(x), label="$X$"); \ + plt.plot(x, Y.pdf(x), label="$Y$"); \ + plt.plot(x, sdr.min_rvs(X, Y).pdf(x), label=r"$\min(X, Y)$"); \ + plt.hist(np.minimum(X.rvs(100_000), Y.rvs(100_000)), bins=101, density=True, histtype="step", label=r"$\min(X, Y)$ empirical"); \ + plt.legend(); \ + plt.xlabel("Random variable"); \ + plt.ylabel("Probability density"); \ + plt.title("Minimum of normal and Rayleigh random variables"); + + Compute the distribution of the minimum of Rayleigh and Rician random variables. + + .. ipython:: python + + X = scipy.stats.rayleigh(scale=2) + Y = scipy.stats.rice(2) + x = np.linspace(0, 8, 1_001) + + @savefig sdr_min_rvs_2.png + plt.figure(); \ + plt.plot(x, X.pdf(x), label="$X$"); \ + plt.plot(x, Y.pdf(x), label="$Y$"); \ + plt.plot(x, sdr.min_rvs(X, Y).pdf(x), label=r"$\min(X, Y)$"); \ + plt.hist(np.minimum(X.rvs(100_000), Y.rvs(100_000)), bins=101, density=True, histtype="step", label=r"$\min(X, Y)$ empirical"); \ + plt.legend(); \ + plt.xlabel("Random variable"); \ + plt.ylabel("Probability density"); \ + plt.title("Minimum of Rayleigh and Rician random variables"); + + Compute the distribution of the minimum of Rician and Chi-squared random variables. + + .. ipython:: python + + X = scipy.stats.rice(3) + Y = scipy.stats.chi2(3) + x = np.linspace(0, 20, 1_001) + + @savefig sdr_min_rvs_2.png + plt.figure(); \ + plt.plot(x, X.pdf(x), label="$X$"); \ + plt.plot(x, Y.pdf(x), label="$Y$"); \ + plt.plot(x, sdr.min_rvs(X, Y).pdf(x), label=r"$\min(X, Y)$"); \ + plt.hist(np.minimum(X.rvs(100_000), Y.rvs(100_000)), bins=101, density=True, histtype="step", label=r"$\min(X, Y)$ empirical"); \ + plt.legend(); \ + plt.xlim(0, 15); \ + plt.xlabel("Random variable"); \ + plt.ylabel("Probability density"); \ + plt.title("Minimum of Rician and Chi-squared random variables"); + + Group: + independent-rvs + """ + verify_scalar(p, float=True, exclusive_min=0, inclusive_max=0.1) + + # Determine the x axis of each distribution such that the probability of exceeding the x axis, on either side, + # is p. + x_min, x_max = _x_range(X, p) + y_min, y_max = _x_range(Y, p) + dx = (x_max - x_min) / 1_000 + dy = (y_max - y_min) / 1_000 + dz = np.min([dx, dy]) # Use the smaller delta x -- must use the same dx for both distributions + z = np.arange(np.min([x_min, y_min]), np.max([x_max, y_max]), dz) + + # Compute the PDF of each distribution + f_X = X.pdf(z) + f_Y = Y.pdf(z) + + # Compute the CDF of each distribution + F_X = X.cdf(z) + F_Y = Y.cdf(z) + + # The PDF of the min of X and Y + f_Z = f_X * (1 - F_Y) + f_Y * (1 - F_X) + + # Adjust the histograms bins to be on either side of each point. So there is one extra point added. + z = np.append(z, z[-1] + dz) + z -= dz / 2 + + return scipy.stats.rv_histogram((f_Z, z)) + + @export def max_rvs( X: scipy.stats.rv_continuous | scipy.stats.rv_histogram, From ae14666fdb22415bd3cc8b1f1a3373a33378e165 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Wed, 18 Dec 2024 06:58:31 -0500 Subject: [PATCH 07/16] Add unit tests for #442 --- tests/probability/test_min_rvs.py | 60 +++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) create mode 100644 tests/probability/test_min_rvs.py diff --git a/tests/probability/test_min_rvs.py b/tests/probability/test_min_rvs.py new file mode 100644 index 000000000..3072f3c6c --- /dev/null +++ b/tests/probability/test_min_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.min_rvs(X, Y) + + # Empirically compute the distribution + z = np.minimum(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("Sum of two distributions") + plt.show() + + assert np.allclose(Z.pdf(x), hist, atol=np.max(hist) * 0.1) From d00f2f293d6a3b038337900cddae69d473de8c39 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Wed, 18 Dec 2024 07:12:47 -0500 Subject: [PATCH 08/16] Fix equations and image names --- src/sdr/_probability.py | 25 +++++++++---------------- 1 file changed, 9 insertions(+), 16 deletions(-) diff --git a/src/sdr/_probability.py b/src/sdr/_probability.py index 308777832..205941ecc 100644 --- a/src/sdr/_probability.py +++ b/src/sdr/_probability.py @@ -428,20 +428,15 @@ def min_rvs( The CDF of $Z$, denoted $F_Z(z)$, is $F_Z(z) = P(Z \leq z)$. Since $Z = \min(X, Y)$, the event $Z \leq z$ occurs if either $X \leq z$ or $Y \leq z$. Using the complement and independence, - $$ - P(Z \leq z) = 1 - P(Z > z) = 1 - P(X > z \text{ and } Y > z) - F_Z(z) = 1 - P(X > z) \cdot P(Y > z) - F_Z(z) = 1 - (1 - F_X(z)) \cdot (1 - F_Y(z)) - $$ + $$P(Z \leq z) = 1 - P(Z > z) = 1 - P(X > z \text{ and } Y > z)$$ + $$F_Z(z) = 1 - P(X > z) \cdot P(Y > z)$$ + $$F_Z(z) = 1 - (1 - F_X(z)) \cdot (1 - 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) = 1 - (1 - F_X(z)) \cdot (1 - F_Y(z))$ yields - $$ - f_Z(z) = \frac{d}{dz} \big(1 - (1 - F_X(z)) \cdot (1 - F_Y(z))\big) - f_Z(z) = \frac{d}{dz} \big((1 - F_X(z)) \cdot (1 - F_Y(z))\big) - f_Z(z) = f_X(z) \cdot (1 - F_Y(z)) + f_Y(z) \cdot (1 - F_X(z)) . - $$ + $$f_Z(z) = \frac{d}{dz} \big(1 - (1 - F_X(z)) \cdot (1 - F_Y(z))\big)$$ + $$f_Z(z) = f_X(z) \cdot (1 - F_Y(z)) + f_Y(z) \cdot (1 - F_X(z)) .$$ Therefore, the PDF of $Z = \min(X, Y)$ is @@ -497,7 +492,7 @@ def min_rvs( Y = scipy.stats.chi2(3) x = np.linspace(0, 20, 1_001) - @savefig sdr_min_rvs_2.png + @savefig sdr_min_rvs_3.png plt.figure(); \ plt.plot(x, X.pdf(x), label="$X$"); \ plt.plot(x, Y.pdf(x), label="$Y$"); \ @@ -572,10 +567,8 @@ def max_rvs( 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) . - $$ + $$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 @@ -631,7 +624,7 @@ def max_rvs( Y = scipy.stats.chi2(3) x = np.linspace(0, 20, 1_001) - @savefig sdr_max_rvs_2.png + @savefig sdr_max_rvs_3.png plt.figure(); \ plt.plot(x, X.pdf(x), label="$X$"); \ plt.plot(x, Y.pdf(x), label="$Y$"); \ From f298eb4db6492f44e5fc72740c58ad92629aa1fd Mon Sep 17 00:00:00 2001 From: mhostetter Date: Thu, 16 Jan 2025 17:19:35 -0500 Subject: [PATCH 09/16] Remove accidental prints Fixes #448 --- src/sdr/_modulation/_cpm.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sdr/_modulation/_cpm.py b/src/sdr/_modulation/_cpm.py index 4d6e5683b..3a5f86170 100644 --- a/src/sdr/_modulation/_cpm.py +++ b/src/sdr/_modulation/_cpm.py @@ -159,10 +159,10 @@ def modulate(self, s: npt.ArrayLike) -> npt.NDArray[np.complex128]: def _modulate(self, s: npt.NDArray[np.int_]) -> npt.NDArray[np.complex128]: s = self._symbol_labels[s] # Relabeled decimal symbols freq = self.index * (2 * s - (self.order - 1)) # Instantaneous frequency - print(freq) + # print(freq) # return f freq_ps = self._tx_pulse_shape(freq) # Pulse-shaped instantaneous frequency - print(freq_ps) + # print(freq_ps) # phase_ps = np.cumsum(freq_ps) # Pulse-shaped instantaneous phase # phase_ps = np.insert(phase_ps, 0, 0) # Start with phase 0 # phase_ps = phase_ps[:-1] # Trim last phase From e4cdf437502335fccb3620c7e11ffd8ee95b3230 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Mon, 20 Jan 2025 23:29:22 -0500 Subject: [PATCH 10/16] Fix bug in `sdr.plot.dft()` Fixes #429 --- src/sdr/plot/_freq_domain.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/sdr/plot/_freq_domain.py b/src/sdr/plot/_freq_domain.py index 7d3faa93f..0706b1037 100644 --- a/src/sdr/plot/_freq_domain.py +++ b/src/sdr/plot/_freq_domain.py @@ -128,8 +128,10 @@ def dft( if ax is None: ax = plt.gca() - if window is not None: - x *= scipy.signal.windows.get_window(window, x.size) + if window is None: + w = np.ones(x.size) + else: + w = scipy.signal.windows.get_window(window, x.size) if size is None: if oversample is None: @@ -140,7 +142,7 @@ def dft( if fast: size = scipy.fft.next_fast_len(size) - X = np.fft.fft(x, size) + X = np.fft.fft(x * w, size) if x_axis == "freq": f = np.fft.fftfreq(size, 1 / sample_rate) elif x_axis == "bin": From 817fe722b0475efe98676e551d41ae6889628787 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Tue, 21 Jan 2025 18:50:06 -0500 Subject: [PATCH 11/16] Bump `galois` version to 0.4.4 --- pyproject.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 3688977e8..b9af9098d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -37,6 +37,7 @@ classifiers = [ "Programming Language :: Python :: 3.10", "Programming Language :: Python :: 3.11", "Programming Language :: Python :: 3.12", + "Programming Language :: Python :: 3.13", "Topic :: Scientific/Engineering :: Mathematics", "Topic :: Security :: Cryptography", "Topic :: Software Development :: Libraries :: Python Modules", @@ -48,7 +49,7 @@ dependencies = [ "numba", # Use galois's version limitation "scipy", "matplotlib", - "galois == 0.4.3", # Exact match required because of use of internals + "galois == 0.4.4", # Exact match required because of use of internals "typing_extensions >= 4.0.0", # v4.0.0 is needed for use of Self (Python 3.11+) and Literal (Python 3.8+) ] dynamic = ["version"] From ed50acf3719d9ea394a76d2f2f5c76b9401e4e2b Mon Sep 17 00:00:00 2001 From: mhostetter Date: Tue, 21 Jan 2025 18:53:32 -0500 Subject: [PATCH 12/16] Test in Python 3.13 --- .github/workflows/lint.yaml | 6 +++--- .github/workflows/test.yaml | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/.github/workflows/lint.yaml b/.github/workflows/lint.yaml index c7e7a9d2f..449c82fe7 100644 --- a/.github/workflows/lint.yaml +++ b/.github/workflows/lint.yaml @@ -16,7 +16,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: [3.8, 3.9, '3.10', 3.11, 3.12] + python-version: [3.8, 3.9, '3.10', 3.11, 3.12, 3.13] steps: - uses: actions/checkout@v3 @@ -43,10 +43,10 @@ jobs: steps: - uses: actions/checkout@v3 - - name: Set up Python 3.10 + - name: Set up Python 3.13 uses: actions/setup-python@v4 with: - python-version: '3.10' + python-version: 3.13 - name: Upgrade pip run: python3 -m pip install --upgrade pip diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index 27f6c07c4..3509e0e4d 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -35,10 +35,10 @@ jobs: steps: - uses: actions/checkout@v3 - - name: Set up Python 3.12 + - name: Set up Python 3.13 uses: actions/setup-python@v4 with: - python-version: 3.12 + python-version: 3.13 - name: Upgrade pip run: python3 -m pip install --upgrade pip @@ -140,7 +140,7 @@ jobs: strategy: matrix: os: [ubuntu-latest, macos-latest, windows-latest] - python-version: [3.8, 3.9, '3.10', 3.11, 3.12] + python-version: [3.8, 3.9, '3.10', 3.11, 3.12, 3.13] runs-on: ${{ matrix.os }} steps: - uses: actions/checkout@v3 From 78b2b47f9d582550db9e79fe0eeac5c373c1b031 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Tue, 21 Jan 2025 18:55:39 -0500 Subject: [PATCH 13/16] Upgrade Ruff to v0.9.2 --- .pre-commit-config.yaml | 2 +- requirements-dev.txt | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index c9150b9b8..1292bc051 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -9,7 +9,7 @@ repos: - id: end-of-file-fixer - id: trailing-whitespace - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.5.1 + rev: v0.9.2 hooks: - id: ruff - id: ruff-format diff --git a/requirements-dev.txt b/requirements-dev.txt index 8f1ac702b..17aa9adfa 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -1,4 +1,4 @@ -ruff == 0.5.1 +ruff == 0.9.2 pre-commit pytest pytest-cov[toml] From 900f61aa1e60307a10747d409496ce84426a9cc5 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Tue, 21 Jan 2025 18:55:49 -0500 Subject: [PATCH 14/16] Fix new linter errors --- src/sdr/_data.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/sdr/_data.py b/src/sdr/_data.py index bf0de44a9..8f08cdddb 100644 --- a/src/sdr/_data.py +++ b/src/sdr/_data.py @@ -153,8 +153,7 @@ def hexdump( data = verify_arraylike(data, int=True, ndim=1, inclusive_min=0, exclusive_max=256) verify_scalar(width, int=True, inclusive_min=1, inclusive_max=16) - if width > data.size: - width = data.size + width = min(width, data.size) string = "" for i in range(0, data.size, width): From e8494c12b95e5f7a7c1d92d2f07e440feb1bd920 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Tue, 21 Jan 2025 18:56:17 -0500 Subject: [PATCH 15/16] Apply new formatting --- src/sdr/_sequence/_correlation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sdr/_sequence/_correlation.py b/src/sdr/_sequence/_correlation.py index fabff417b..da20a4255 100644 --- a/src/sdr/_sequence/_correlation.py +++ b/src/sdr/_sequence/_correlation.py @@ -627,7 +627,7 @@ def _kasami_small_set(degree: int, poly: Poly, index: int) -> npt.NDArray[np.int j = index if not -1 <= j < 2 ** (degree // 2) - 1: - raise ValueError(f"Argument 'index' must be between -1 and {2**(degree//2) - 1}, not {j}.") + raise ValueError(f"Argument 'index' must be between -1 and {2 ** (degree // 2) - 1}, not {j}.") u = m_sequence(degree, poly=poly, index=1, output="decimal") length = u.size @@ -652,7 +652,7 @@ def _kasami_large_set(degree: int, poly: Poly, index: tuple[int, int]) -> npt.ND if not -2 <= i < 2**degree - 1: raise ValueError(f"Argument 'index[0]' must be between -2 and {2**degree - 1}, not {i}.") if not -1 <= j < 2 ** (degree // 2) - 1: - raise ValueError(f"Argument 'index[1]' must be between -1 and {2**(degree//2) - 1}, not {j}.") + raise ValueError(f"Argument 'index[1]' must be between -1 and {2 ** (degree // 2) - 1}, not {j}.") u = m_sequence(degree, poly=poly, index=1, output="decimal") length = u.size From 9a5b649448255a3edd58d7e7858841ae33b7826a Mon Sep 17 00:00:00 2001 From: mhostetter Date: Tue, 21 Jan 2025 18:58:16 -0500 Subject: [PATCH 16/16] Add release notes for v0.0.26 --- docs/release-notes/v0.0.md | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/docs/release-notes/v0.0.md b/docs/release-notes/v0.0.md index 1a5061449..5c1a9f2d7 100644 --- a/docs/release-notes/v0.0.md +++ b/docs/release-notes/v0.0.md @@ -4,6 +4,21 @@ tocdepth: 2 # v0.0 +## v0.0.26 + +*Released January 21, 2025* + +### Changes + +- Pinned `galois` version to 0.4.4. +- Added support for Python 3.13. +- Added support for NumPy 2.1. +- Added support for Numba 0.61. + +### Contributors + +- Matt Hostetter ([@mhostetter](https://github.com/mhostetter)) + ## v0.0.25 *Released December 16, 2024*