From 31fef21e9547ab9f88b8a9ad1b48793d6e0e3ee2 Mon Sep 17 00:00:00 2001 From: Colm Talbot Date: Wed, 10 Jun 2026 10:30:50 -0400 Subject: [PATCH 1/4] FEAT: add skew t distribution --- gwpopulation/utils.py | 36 ++++++++++++++++++++++++++++++++++++ test/utils_test.py | 22 +++++++++++++++++++++- 2 files changed, 57 insertions(+), 1 deletion(-) diff --git a/gwpopulation/utils.py b/gwpopulation/utils.py index b8a909af..21e77fc7 100644 --- a/gwpopulation/utils.py +++ b/gwpopulation/utils.py @@ -183,6 +183,42 @@ def logsubexp(log_p, log_q): return xp.nan_to_num(xp.exp(log_pdf)) * (xx >= low) * (xx <= high) + +def skewt(xx, aa, bb, loc=0, scale=1): + r""" + Jones and Faddy skew-t distribution (implementation based on :code:`scipy`). + + .. math:: + + z &= \frac{x - \mu}{\sigma} \\ + p(x) &= \frac{1}{\sigma C_{a,b}} \left(1 + \frac{z}{\sqrt{a + b + z^2}}\right)^{a + \frac{1}{2}} \left(1 - \frac{z}{\sqrt{a + b + z^2}}\right)^{b + \frac{1}{2}} \\ + C_{a,b} &= {2^{a + b - 1} B(a, b) \sqrt{a + b}} + + Parameters + ---------- + xx: float, array-like + The abscissa values (:math:`x`) + aa: float + The first shape parameter (:math:`a`) + bb: float + The second shape parameter (:math:`b`) + loc: float + The location parameter (:math:`\mu`) + scale: float + The scale parameter (:math:`\sigma`) + + Returns + ------- + prob: float, array-like + The distribution evaluated at `xx` + """ + zz = (xx - loc) / scale + c = 2 ** (aa + bb - 1) * scs.beta(aa, bb) * (aa + bb)**0.5 + d1 = (1 + zz / (aa + bb + zz ** 2)**0.5) ** (aa + 0.5) + d2 = (1 - zz / (aa + bb + zz ** 2)**0.5) ** (bb + 0.5) + return d1 * d2 / c / scale + + def unnormalized_2d_gaussian(xx, yy, mu_x, mu_y, sigma_x, sigma_y, covariance): r""" Compute the probability distribution for a correlated 2-dimensional Gaussian diff --git a/test/utils_test.py b/test/utils_test.py index 56a7c66c..98fb5585 100644 --- a/test/utils_test.py +++ b/test/utils_test.py @@ -127,7 +127,27 @@ def test_truncnorm_matches_scipy(backend): assert max(abs(gwpop_vals - scipy_vals)) < 1e-3 -def test_matches_scipy(backend): +def test_skewt_matches_scipy(backend): + from scipy.stats import fj_skew_t + + gwpopulation.set_backend(backend) + xp = gwpopulation.utils.xp + xx = xp.linspace(-2, 2, 1000) + for ii in range(N_TEST): + mu = np.random.uniform(-10, 10) + sigma = np.random.uniform(0, 5) + aa = np.random.uniform(0, 100) + bb = np.random.uniform(0, 100) + gwpop_vals = utils.to_numpy( + utils.skewt(xx, aa=aa, bb=bb, loc=mu, scale=sigma) + ) + scipy_vals = fj_skew_t( + loc=mu, scale=sigma, a=aa, b=bb + ).pdf(utils.to_numpy(xx)) + assert max(abs(gwpop_vals - scipy_vals)) < 1e-3 + + +def test_vonmises_matches_scipy(backend): gwpopulation.set_backend(backend) xp = gwpopulation.utils.xp xx = xp.linspace(0, 2 * np.pi, 1000) From 43db0f442a46ecf2738513d355c9bf1fb5849d95 Mon Sep 17 00:00:00 2001 From: Colm Talbot Date: Wed, 10 Jun 2026 10:48:54 -0400 Subject: [PATCH 2/4] Apply suggestions from code review Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> Co-authored-by: Colm Talbot --- gwpopulation/utils.py | 20 ++++++++++++++------ test/utils_test.py | 4 ++-- 2 files changed, 16 insertions(+), 8 deletions(-) diff --git a/gwpopulation/utils.py b/gwpopulation/utils.py index 21e77fc7..2c6ed030 100644 --- a/gwpopulation/utils.py +++ b/gwpopulation/utils.py @@ -183,7 +183,7 @@ def logsubexp(log_p, log_q): return xp.nan_to_num(xp.exp(log_pdf)) * (xx >= low) * (xx <= high) - +@apply_conditions(dict(aa=(gt, 0), bb=(gt, 0), scale=(gt, 0))) def skewt(xx, aa, bb, loc=0, scale=1): r""" Jones and Faddy skew-t distribution (implementation based on :code:`scipy`). @@ -212,11 +212,19 @@ def skewt(xx, aa, bb, loc=0, scale=1): prob: float, array-like The distribution evaluated at `xx` """ - zz = (xx - loc) / scale - c = 2 ** (aa + bb - 1) * scs.beta(aa, bb) * (aa + bb)**0.5 - d1 = (1 + zz / (aa + bb + zz ** 2)**0.5) ** (aa + 0.5) - d2 = (1 - zz / (aa + bb + zz ** 2)**0.5) ** (bb + 0.5) - return d1 * d2 / c / scale + denom = xp.sqrt(aa + bb + zz**2) + log_c = ( + (aa + bb - 1) * np.log(2) + + scs.betaln(aa, bb) + + xp.log(aa + bb) / 2 + + xp.log(scale) + ) + log_pdf = ( + (aa + 0.5) * xp.log1p(zz / denom) + + (bb + 0.5) * xp.log1p(-zz / denom) + - log_c + ) + return xp.nan_to_num(xp.exp(log_pdf)) def unnormalized_2d_gaussian(xx, yy, mu_x, mu_y, sigma_x, sigma_y, covariance): diff --git a/test/utils_test.py b/test/utils_test.py index 98fb5585..9546207b 100644 --- a/test/utils_test.py +++ b/test/utils_test.py @@ -128,7 +128,7 @@ def test_truncnorm_matches_scipy(backend): def test_skewt_matches_scipy(backend): - from scipy.stats import fj_skew_t + from scipy.stats import jf_skew_t gwpopulation.set_backend(backend) xp = gwpopulation.utils.xp @@ -141,7 +141,7 @@ def test_skewt_matches_scipy(backend): gwpop_vals = utils.to_numpy( utils.skewt(xx, aa=aa, bb=bb, loc=mu, scale=sigma) ) - scipy_vals = fj_skew_t( + scipy_vals = jf_skew_t( loc=mu, scale=sigma, a=aa, b=bb ).pdf(utils.to_numpy(xx)) assert max(abs(gwpop_vals - scipy_vals)) < 1e-3 From fbdfa8ceae56a62b2bb60f8929e60781a8b33def Mon Sep 17 00:00:00 2001 From: Colm Talbot Date: Wed, 10 Jun 2026 10:53:27 -0400 Subject: [PATCH 3/4] TYPO: add removed line --- gwpopulation/utils.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/gwpopulation/utils.py b/gwpopulation/utils.py index 2c6ed030..16e13a3a 100644 --- a/gwpopulation/utils.py +++ b/gwpopulation/utils.py @@ -191,7 +191,9 @@ def skewt(xx, aa, bb, loc=0, scale=1): .. math:: z &= \frac{x - \mu}{\sigma} \\ - p(x) &= \frac{1}{\sigma C_{a,b}} \left(1 + \frac{z}{\sqrt{a + b + z^2}}\right)^{a + \frac{1}{2}} \left(1 - \frac{z}{\sqrt{a + b + z^2}}\right)^{b + \frac{1}{2}} \\ + p(x) &= \frac{1}{\sigma C_{a,b}} + \left(1 + \frac{z}{\sqrt{a + b + z^2}}\right)^{a + \frac{1}{2}} + \left(1 - \frac{z}{\sqrt{a + b + z^2}}\right)^{b + \frac{1}{2}} \\ C_{a,b} &= {2^{a + b - 1} B(a, b) \sqrt{a + b}} Parameters @@ -212,6 +214,7 @@ def skewt(xx, aa, bb, loc=0, scale=1): prob: float, array-like The distribution evaluated at `xx` """ + zz = (xx - loc) / scale denom = xp.sqrt(aa + bb + zz**2) log_c = ( (aa + bb - 1) * np.log(2) From db99fab7c0d1aafd4fa9a84514ed379e36d6a74d Mon Sep 17 00:00:00 2001 From: Colm Talbot Date: Wed, 10 Jun 2026 11:03:39 -0400 Subject: [PATCH 4/4] FMT: fix formatting --- gwpopulation/utils.py | 4 +--- test/utils_test.py | 8 ++------ 2 files changed, 3 insertions(+), 9 deletions(-) diff --git a/gwpopulation/utils.py b/gwpopulation/utils.py index 16e13a3a..8d1b0b53 100644 --- a/gwpopulation/utils.py +++ b/gwpopulation/utils.py @@ -223,9 +223,7 @@ def skewt(xx, aa, bb, loc=0, scale=1): + xp.log(scale) ) log_pdf = ( - (aa + 0.5) * xp.log1p(zz / denom) - + (bb + 0.5) * xp.log1p(-zz / denom) - - log_c + (aa + 0.5) * xp.log1p(zz / denom) + (bb + 0.5) * xp.log1p(-zz / denom) - log_c ) return xp.nan_to_num(xp.exp(log_pdf)) diff --git a/test/utils_test.py b/test/utils_test.py index 9546207b..7e4e8d75 100644 --- a/test/utils_test.py +++ b/test/utils_test.py @@ -138,12 +138,8 @@ def test_skewt_matches_scipy(backend): sigma = np.random.uniform(0, 5) aa = np.random.uniform(0, 100) bb = np.random.uniform(0, 100) - gwpop_vals = utils.to_numpy( - utils.skewt(xx, aa=aa, bb=bb, loc=mu, scale=sigma) - ) - scipy_vals = jf_skew_t( - loc=mu, scale=sigma, a=aa, b=bb - ).pdf(utils.to_numpy(xx)) + gwpop_vals = utils.to_numpy(utils.skewt(xx, aa=aa, bb=bb, loc=mu, scale=sigma)) + scipy_vals = jf_skew_t(loc=mu, scale=sigma, a=aa, b=bb).pdf(utils.to_numpy(xx)) assert max(abs(gwpop_vals - scipy_vals)) < 1e-3