From 0dcd76c55c09606337e89e203d28a921f2c341cc Mon Sep 17 00:00:00 2001 From: Matt Haberland Date: Sun, 9 Feb 2025 12:16:31 -0800 Subject: [PATCH] MAINT: `stats.invgauss`: return correct result when `mu=inf` (#22496) --- scipy/stats/_continuous_distns.py | 12 ++++++------ scipy/stats/tests/test_distributions.py | 16 ++++++++++++++++ 2 files changed, 22 insertions(+), 6 deletions(-) diff --git a/scipy/stats/_continuous_distns.py b/scipy/stats/_continuous_distns.py index dfd517b1e0fc..4220680dc38c 100644 --- a/scipy/stats/_continuous_distns.py +++ b/scipy/stats/_continuous_distns.py @@ -5018,10 +5018,10 @@ def _rvs(self, mu, size=None, random_state=None): def _pdf(self, x, mu): # invgauss.pdf(x, mu) = # 1 / sqrt(2*pi*x**3) * exp(-(x-mu)**2/(2*x*mu**2)) - return 1.0/np.sqrt(2*np.pi*x**3.0)*np.exp(-1.0/(2*x)*((x-mu)/mu)**2) + return 1.0/np.sqrt(2*np.pi*x**3.0)*np.exp(-1.0/(2*x)*(x/mu - 1)**2) def _logpdf(self, x, mu): - return -0.5*np.log(2*np.pi) - 1.5*np.log(x) - ((x-mu)/mu)**2/(2*x) + return -0.5*np.log(2*np.pi) - 1.5*np.log(x) - (x/mu - 1)**2/(2*x) # approach adapted from equations in # https://journal.r-project.org/archive/2016-1/giner-smyth.pdf, @@ -5029,14 +5029,14 @@ def _logpdf(self, x, mu): def _logcdf(self, x, mu): fac = 1 / np.sqrt(x) - a = _norm_logcdf(fac * ((x / mu) - 1)) - b = 2 / mu + _norm_logcdf(-fac * ((x / mu) + 1)) + a = _norm_logcdf(fac * (x/mu - 1)) + b = 2 / mu + _norm_logcdf(-fac * (x/mu + 1)) return a + np.log1p(np.exp(b - a)) def _logsf(self, x, mu): fac = 1 / np.sqrt(x) - a = _norm_logsf(fac * ((x / mu) - 1)) - b = 2 / mu + _norm_logcdf(-fac * (x + mu) / mu) + a = _norm_logsf(fac * (x/mu - 1)) + b = 2 / mu + _norm_logcdf(-fac * (x/mu + 1)) return a + np.log1p(-np.exp(b - a)) def _sf(self, x, mu): diff --git a/scipy/stats/tests/test_distributions.py b/scipy/stats/tests/test_distributions.py index 71df924f1d15..32d800323d33 100644 --- a/scipy/stats/tests/test_distributions.py +++ b/scipy/stats/tests/test_distributions.py @@ -3345,6 +3345,22 @@ def test_logcdf_logsf(self): def test_entropy(self, mu, ref): assert_allclose(stats.invgauss.entropy(mu), ref, rtol=5e-14) + def test_mu_inf_gh13666(self): + # invgauss methods should return correct result when mu=inf + # invgauss as mu -> oo is invgamma with shape and scale 0.5; + # see gh-13666 and gh-22496 + dist = stats.invgauss(mu=np.inf) + dist0 = stats.invgamma(0.5, scale=0.5) + x, p = 1., 0.5 + assert_allclose(dist.logpdf(x), dist0.logpdf(x)) + assert_allclose(dist.pdf(x), dist0.pdf(x)) + assert_allclose(dist.logcdf(x), dist0.logcdf(x)) + assert_allclose(dist.cdf(x), dist0.cdf(x)) + assert_allclose(dist.logsf(x), dist0.logsf(x)) + assert_allclose(dist.sf(x), dist0.sf(x)) + assert_allclose(dist.ppf(p), dist0.ppf(p)) + assert_allclose(dist.isf(p), dist0.isf(p)) + class TestLandau: @pytest.mark.parametrize('name', ['pdf', 'cdf', 'sf', 'ppf', 'isf'])