From f69ac4da3063bff083ddaa89570295d6fbc7db1f Mon Sep 17 00:00:00 2001 From: Thomas Wiecki Date: Sun, 5 Mar 2017 02:39:41 +0100 Subject: [PATCH 1/4] Fixes to Rice test. --- pymc3/distributions/continuous.py | 75 ++++++++++++++++++++++++++----- pymc3/tests/test_distributions.py | 4 ++ 2 files changed, 68 insertions(+), 11 deletions(-) diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index f0e63e5e68..ce9b4719b5 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -26,7 +26,8 @@ 'Laplace', 'StudentT', 'Cauchy', 'HalfCauchy', 'Gamma', 'Weibull', 'HalfStudentT', 'Lognormal', 'ChiSquared', 'HalfNormal', 'Wald', 'Pareto', 'InverseGamma', 'ExGaussian', 'VonMises', 'SkewNormal', - 'Triangular', 'Gumbel', 'Logistic', 'LogitNormal', 'Interpolated'] + 'Triangular', 'Gumbel', 'Logistic', 'LogitNormal', 'Interpolated', + 'Rice'] class PositiveContinuous(Continuous): @@ -68,7 +69,7 @@ def assert_negative_support(var, label, distname, value=-1e-6): def get_tau_sd(tau=None, sd=None): """ - Find precision and standard deviation. The link between the two + Find precision and standard deviation. The link between the two parameterizations is given by the inverse relationship: .. math:: @@ -275,14 +276,14 @@ class Normal(Continuous): Standard deviation (sd > 0) (only required if tau is not specified). tau : float Precision (tau > 0) (only required if sd is not specified). - + Examples -------- .. code-block:: python with pm.Model(): x = pm.Normal('x', mu=0, sd=10) - + with pm.Model(): x = pm.Normal('x', mu=0, tau=1/23) """ @@ -887,7 +888,7 @@ class Lognormal(PositiveContinuous): Standard deviation. (sd > 0). (only required if tau is not specified). tau : float Scale parameter (tau > 0). (only required if sd is not specified). - + Example ------- .. code-block:: python @@ -994,14 +995,14 @@ class StudentT(Continuous): Standard deviation (sd > 0) (only required if lam is not specified) lam : float Precision (lam > 0) (only required if sd is not specified) - + Examples -------- .. code-block:: python with pm.Model(): x = pm.StudentT('x', nu=15, mu=0, sd=10) - + with pm.Model(): x = pm.StudentT('x', nu=15, mu=0, lam=1/23) """ @@ -1601,8 +1602,8 @@ def __init__(self, alpha, beta, *args, **kwargs): self.median = beta * tt.exp(gammaln(tt.log(2)))**(1. / alpha) self.variance = (beta**2) * \ tt.exp(gammaln(1 + 2. / alpha - self.mean**2)) - self.mode = tt.switch(alpha >= 1, - beta * ((alpha - 1)/alpha) ** (1 / alpha), + self.mode = tt.switch(alpha >= 1, + beta * ((alpha - 1)/alpha) ** (1 / alpha), 0) # Reference: https://en.wikipedia.org/wiki/Weibull_distribution assert_negative_support(alpha, 'alpha', 'Weibull') @@ -1682,7 +1683,7 @@ class HalfStudentT(PositiveContinuous): lam : float Scale parameter (lam > 0). Converges to the precision as nu increases. (only required if sd is not specified) - + Examples -------- .. code-block:: python @@ -1690,7 +1691,7 @@ class HalfStudentT(PositiveContinuous): # Only pass in one of lam or sd, but not both. with pm.Model(): x = pm.HalfStudentT('x', sd=10, nu=10) - + with pm.Model(): x = pm.HalfStudentT('x', lam=4, nu=10) """ @@ -2198,6 +2199,58 @@ def _repr_latex_(self, name=None, dist=None): get_variable_name(beta)) +class Rice(Continuous): + R""" + Rice distribution. + + .. math:: + + f(x\mid \nu ,\sigma )= + {\frac {x}{\sigma ^{2}}}\exp + \left({\frac {-(x^{2}+\nu ^{2})}{2\sigma ^{2}}}\right)I_{0}\left({\frac {x\nu }{\sigma ^{2}}}\right), + + ======== ============================================================== + Support :math:`x \in (0, +\infinity)` + Mean :math:`\sigma {\sqrt {\pi /2}}\,\,L_{{1/2}}(-\nu ^{2}/2\sigma ^{2})` + Variance :math:`2\sigma ^{2}+\nu ^{2}-{\frac {\pi \sigma ^{2}}{2}}L_{{1/2}}^{2} + \left({\frac {-\nu ^{2}}{2\sigma ^{2}}}\right)` + ======== ============================================================== + + + Parameters + ---------- + nu : float + shape parameter. + sd : float + standard deviation. + + """ + def __init__(self, nu=None, sd=None, *args, **kwargs): + super(Rice, self).__init__(*args, **kwargs) + self.nu = nu = tt.as_tensor_variable(nu) + self.sd = sd = tt.as_tensor_variable(sd) + self.mean = sd * np.sqrt(np.pi / 2) * tt.exp((-nu**2 / (2 * sd**2)) / 2) * ((1 - (-nu**2 / (2 * sd**2))) + * i0(-(-nu**2 / (2 * sd**2)) / 2) - (-nu**2 / (2 * sd**2)) * i1(-(-nu**2 / (2 * sd**2)) / 2)) + self.variance = 2 * sd**2 + nu**2 - (np.pi * sd**2 / 2) * (tt.exp((-nu**2 / (2 * sd**2)) / 2) * ((1 - (-nu**2 / ( + 2 * sd**2))) * i0(-(-nu**2 / (2 * sd**2)) / 2) - (-nu**2 / (2 * sd**2)) * i1(-(-nu**2 / (2 * sd**2)) / 2)))**2 + + + def random(self, point=None, size=None, repeat=None): + nu, sd = draw_values([self.nu, self.sd], + point=point) + return generate_samples(stats.rice.rvs, b=nu, scale=sd, loc=0, + dist_shape=self.shape, size=size) + + def logp(self, value): + nu = self.nu + sd = self.sd + return bound(tt.log(value / (sd**2) * tt.exp(-(value**2 + nu**2) / (2 * sd**2)) * i0(value * nu / (sd**2))), + sd >= 0, + nu >= 0, + value > 0, + ) + + class Logistic(Continuous): R""" Logistic log-likelihood. diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index ccee66a69a..d8702d3ffe 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -1056,6 +1056,10 @@ def test_multidimensional_beta_construction(self): with Model(): Beta('beta', alpha=1., beta=1., shape=(10, 20)) + def test_rice(self): + self.pymc3_matches_scipy(Rice, Rplus, {'nu': Rplus, 'sd': Rplus}, + lambda value, nu, sd: sp.rice.logpdf(value / sd, loc=nu, scale=sd)) + @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32") def test_interpolated(self): for mu in R.vals: From 4512a6b230e2c59d474aeed422ed3bd4f2948a83 Mon Sep 17 00:00:00 2001 From: MBrouns Date: Wed, 2 May 2018 20:58:47 +0200 Subject: [PATCH 2/4] make test case for Rice distribution pass --- pymc3/distributions/__init__.py | 1 + pymc3/distributions/continuous.py | 9 +++++---- pymc3/distributions/dist_math.py | 20 ++++++++++++++++++++ pymc3/tests/test_distributions.py | 7 ++++--- 4 files changed, 30 insertions(+), 7 deletions(-) diff --git a/pymc3/distributions/__init__.py b/pymc3/distributions/__init__.py index 3f1c2b6a29..6d7abce67f 100644 --- a/pymc3/distributions/__init__.py +++ b/pymc3/distributions/__init__.py @@ -28,6 +28,7 @@ from .continuous import Logistic from .continuous import LogitNormal from .continuous import Interpolated +from .continuous import Rice from .discrete import Binomial from .discrete import BetaBinomial diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index ce9b4719b5..845715a63c 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -19,7 +19,7 @@ from pymc3.util import get_variable_name from .special import log_i0 from ..math import invlogit, logit -from .dist_math import bound, logpow, gammaln, betaln, std_cdf, alltrue_elemwise, SplineWrapper +from .dist_math import bound, logpow, gammaln, betaln, std_cdf, alltrue_elemwise, SplineWrapper, i0, i1 from .distribution import Continuous, draw_values, generate_samples __all__ = ['Uniform', 'Flat', 'HalfFlat', 'Normal', 'Beta', 'Exponential', @@ -2237,18 +2237,19 @@ def __init__(self, nu=None, sd=None, *args, **kwargs): def random(self, point=None, size=None, repeat=None): nu, sd = draw_values([self.nu, self.sd], - point=point) + point=point) return generate_samples(stats.rice.rvs, b=nu, scale=sd, loc=0, dist_shape=self.shape, size=size) def logp(self, value): nu = self.nu sd = self.sd - return bound(tt.log(value / (sd**2) * tt.exp(-(value**2 + nu**2) / (2 * sd**2)) * i0(value * nu / (sd**2))), + x = value / sd + return bound(tt.log(x * tt.exp((-(x ** 2 + nu ** 2)) / 2)) + log_i0(x * nu) - tt.log(sd), sd >= 0, nu >= 0, value > 0, - ) + ) class Logistic(Continuous): diff --git a/pymc3/distributions/dist_math.py b/pymc3/distributions/dist_math.py index e6e6daa1cf..1184ec1972 100644 --- a/pymc3/distributions/dist_math.py +++ b/pymc3/distributions/dist_math.py @@ -461,3 +461,23 @@ def infer_shape(self, node, shapes): def eigh(a, UPLO='L'): """A copy, remove with Eigh and EighGrad when possible""" return Eigh(UPLO)(a) + + +def i0(x): + """ + Calculates the 0 order modified Bessel function of the first kind"" + """ + return tt.switch(tt.lt(x, 5), 1 + x**2 / 4 + x**4 / 64 + x**6 / 2304 + x**8 / 147456 + + x**10 / 14745600 + x**12 / 2123366400, + np.e**x / (2 * np.pi * x)**0.5 * (1 + 1 / (8 * x) + 9 / (128 * x**2) + 225 / (3072 * x**3) + + 11025 / (98304 * x**4))) + + +def i1(x): + """ + Calculates the 1 order modified Bessel function of the first kind"" + """ + return tt.switch(tt.lt(x, 5), x / 2 + x**3 / 16 + x**5 / 384 + x**7 / 18432 + + x**9 / 1474560 + x**11 / 176947200 + x**13 / 29727129600, + np.e**x / (2 * np.pi * x)**0.5 * (1 - 3 / (8 * x) + 15 / (128 * x**2) + 315 / (3072 * x**3) + + 14175 / (98304 * x**4))) diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index d8702d3ffe..152028d379 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -16,7 +16,7 @@ Flat, LKJCorr, Wald, ChiSquared, HalfNormal, DiscreteUniform, Bound, Uniform, Triangular, Binomial, SkewNormal, DiscreteWeibull, Gumbel, Logistic, OrderedLogistic, LogitNormal, Interpolated, - ZeroInflatedBinomial, HalfFlat, AR1, KroneckerNormal) + ZeroInflatedBinomial, HalfFlat, AR1, KroneckerNormal, Rice) from ..distributions import continuous from pymc3.theanof import floatX @@ -1057,8 +1057,9 @@ def test_multidimensional_beta_construction(self): Beta('beta', alpha=1., beta=1., shape=(10, 20)) def test_rice(self): - self.pymc3_matches_scipy(Rice, Rplus, {'nu': Rplus, 'sd': Rplus}, - lambda value, nu, sd: sp.rice.logpdf(value / sd, loc=nu, scale=sd)) + self.pymc3_matches_scipy(Rice, Rplusbig, {'nu': Rplusbig, 'sd': Rplusbig}, + lambda value, nu, sd: sp.rice.logpdf(value, b=nu, loc=0, scale=sd), + decimal=select_by_precision(float64=3, float32=1)) @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32") def test_interpolated(self): From 8f87ba203bb60b53b55ba103c88b9d57a87a37d3 Mon Sep 17 00:00:00 2001 From: MBrouns Date: Thu, 17 May 2018 19:08:42 +0200 Subject: [PATCH 3/4] Improve numerical stability of rice distribution function for large values of nu and value Use ((x-b)**2)/2 + xb instead of (x**2 + b**2)/2 in the pdf for the rice distribution and include the np.exp(-xb) in the i0e to match the scipy implementation --- pymc3/distributions/continuous.py | 8 +++++--- pymc3/distributions/dist_math.py | 22 ++++++++-------------- pymc3/tests/test_distributions.py | 4 ++-- 3 files changed, 15 insertions(+), 19 deletions(-) diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index 845715a63c..9272ad42d3 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -14,12 +14,14 @@ from scipy.interpolate import InterpolatedUnivariateSpline import warnings +from theano.scalar import i1, i0 + from pymc3.theanof import floatX from . import transforms from pymc3.util import get_variable_name from .special import log_i0 from ..math import invlogit, logit -from .dist_math import bound, logpow, gammaln, betaln, std_cdf, alltrue_elemwise, SplineWrapper, i0, i1 +from .dist_math import bound, logpow, gammaln, betaln, std_cdf, alltrue_elemwise, SplineWrapper, i0e from .distribution import Continuous, draw_values, generate_samples __all__ = ['Uniform', 'Flat', 'HalfFlat', 'Normal', 'Beta', 'Exponential', @@ -2225,6 +2227,7 @@ class Rice(Continuous): standard deviation. """ + def __init__(self, nu=None, sd=None, *args, **kwargs): super(Rice, self).__init__(*args, **kwargs) self.nu = nu = tt.as_tensor_variable(nu) @@ -2234,7 +2237,6 @@ def __init__(self, nu=None, sd=None, *args, **kwargs): self.variance = 2 * sd**2 + nu**2 - (np.pi * sd**2 / 2) * (tt.exp((-nu**2 / (2 * sd**2)) / 2) * ((1 - (-nu**2 / ( 2 * sd**2))) * i0(-(-nu**2 / (2 * sd**2)) / 2) - (-nu**2 / (2 * sd**2)) * i1(-(-nu**2 / (2 * sd**2)) / 2)))**2 - def random(self, point=None, size=None, repeat=None): nu, sd = draw_values([self.nu, self.sd], point=point) @@ -2245,7 +2247,7 @@ def logp(self, value): nu = self.nu sd = self.sd x = value / sd - return bound(tt.log(x * tt.exp((-(x ** 2 + nu ** 2)) / 2)) + log_i0(x * nu) - tt.log(sd), + return bound(tt.log(x * tt.exp((-(x - nu) * (x - nu)) / 2) * i0e(x * nu) / sd), sd >= 0, nu >= 0, value > 0, diff --git a/pymc3/distributions/dist_math.py b/pymc3/distributions/dist_math.py index 1184ec1972..f6725953cc 100644 --- a/pymc3/distributions/dist_math.py +++ b/pymc3/distributions/dist_math.py @@ -9,6 +9,7 @@ import scipy.linalg import theano.tensor as tt import theano +from theano.scalar import UnaryScalarOp, upgrade_to_float from .special import gammaln from pymc3.theanof import floatX @@ -463,21 +464,14 @@ def eigh(a, UPLO='L'): return Eigh(UPLO)(a) -def i0(x): +class I0e(UnaryScalarOp): """ - Calculates the 0 order modified Bessel function of the first kind"" + Modified Bessel function of the first kind of order 0, exponentially scaled. """ - return tt.switch(tt.lt(x, 5), 1 + x**2 / 4 + x**4 / 64 + x**6 / 2304 + x**8 / 147456 - + x**10 / 14745600 + x**12 / 2123366400, - np.e**x / (2 * np.pi * x)**0.5 * (1 + 1 / (8 * x) + 9 / (128 * x**2) + 225 / (3072 * x**3) - + 11025 / (98304 * x**4))) + nfunc_spec = ('scipy.special.i0e', 1, 1) + def impl(self, x): + return scipy.special.i0e(x) -def i1(x): - """ - Calculates the 1 order modified Bessel function of the first kind"" - """ - return tt.switch(tt.lt(x, 5), x / 2 + x**3 / 16 + x**5 / 384 + x**7 / 18432 + - x**9 / 1474560 + x**11 / 176947200 + x**13 / 29727129600, - np.e**x / (2 * np.pi * x)**0.5 * (1 - 3 / (8 * x) + 15 / (128 * x**2) + 315 / (3072 * x**3) - + 14175 / (98304 * x**4))) + +i0e = I0e(upgrade_to_float, name='i0e') diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index 152028d379..cb8f8a12d7 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -1057,9 +1057,9 @@ def test_multidimensional_beta_construction(self): Beta('beta', alpha=1., beta=1., shape=(10, 20)) def test_rice(self): - self.pymc3_matches_scipy(Rice, Rplusbig, {'nu': Rplusbig, 'sd': Rplusbig}, + self.pymc3_matches_scipy(Rice, Rplus, {'nu': Rplus, 'sd': Rplus}, lambda value, nu, sd: sp.rice.logpdf(value, b=nu, loc=0, scale=sd), - decimal=select_by_precision(float64=3, float32=1)) + decimal=select_by_precision(float64=6, float32=1)) @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32") def test_interpolated(self): From ed1fac9f82c47de5435f260253384589bf799a98 Mon Sep 17 00:00:00 2001 From: MBrouns Date: Sat, 19 May 2018 15:33:37 +0200 Subject: [PATCH 4/4] Change test domain for sd of rice distribution to pass tests with float32 --- pymc3/tests/test_distributions.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index cb8f8a12d7..8ef4e0afe0 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -1057,9 +1057,8 @@ def test_multidimensional_beta_construction(self): Beta('beta', alpha=1., beta=1., shape=(10, 20)) def test_rice(self): - self.pymc3_matches_scipy(Rice, Rplus, {'nu': Rplus, 'sd': Rplus}, - lambda value, nu, sd: sp.rice.logpdf(value, b=nu, loc=0, scale=sd), - decimal=select_by_precision(float64=6, float32=1)) + self.pymc3_matches_scipy(Rice, Rplus, {'nu': Rplus, 'sd': Rplusbig}, + lambda value, nu, sd: sp.rice.logpdf(value, b=nu, loc=0, scale=sd)) @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32") def test_interpolated(self):