From e4e8734d8d2771ef39edbb2ecc8bb33647842b70 Mon Sep 17 00:00:00 2001 From: "M. Domenzain" Date: Mon, 6 Mar 2017 19:04:10 +0100 Subject: [PATCH 1/7] Implement log CDF for Uniform distribution --- pymc3/distributions/continuous.py | 12 ++++++++++++ pymc3/tests/test_distributions.py | 2 ++ 2 files changed, 14 insertions(+) diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index 746d4aa240..b639d2142d 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -152,6 +152,18 @@ def logp(self, value): return bound(-tt.log(upper - lower), value >= lower, value <= upper) + def logcdf(self, value): + return tt.switch( + tt.or_(tt.lt(value, self.lower), tt.gt(value, self.upper)), + -np.inf, + tt.switch( + tt.eq(value, self.upper), + 0, + tt.log((value - self.lower)) - + tt.log((self.upper - self.lower)) + ) + ) + class Flat(Continuous): """ diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index 9449bb8a78..f06402b4ca 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -366,6 +366,8 @@ def test_uniform(self): self.pymc3_matches_scipy( Uniform, Runif, {'lower': -Rplusunif, 'upper': Rplusunif}, lambda value, lower, upper: sp.uniform.logpdf(value, lower, upper - lower)) + self.check_logcdf(Uniform, Runif, {'lower': -Rplusunif, 'upper': Rplusunif}, + lambda value, lower, upper: sp.uniform.logcdf(value, lower, upper - lower)) def test_triangular(self): self.pymc3_matches_scipy( From 4fb89a08cb87cc1decbc6e1bb937f40773255a55 Mon Sep 17 00:00:00 2001 From: "M. Domenzain" Date: Mon, 6 Mar 2017 19:08:22 +0100 Subject: [PATCH 2/7] Implement log CDF for HalfNormal distribution --- pymc3/distributions/continuous.py | 11 ++++++++++- pymc3/tests/test_distributions.py | 4 +++- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index b639d2142d..1d8153d6f7 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -255,7 +255,7 @@ def logcdf(self, value): tt.sqr(tt.abs_(z)) / 2, tt.log1p(-tt.erfc(z / tt.sqrt(2.)) / 2.) ) - + class HalfNormal(PositiveContinuous): R""" @@ -307,6 +307,15 @@ def logp(self, value): value >= 0, tau > 0, sd > 0) + def logcdf(self, value): + sd = self.sd + z = zvalue(value, mu=0, sd=sd) + return tt.switch( + tt.lt(z, -1.0), + tt.log(tt.erfcx(-z / tt.sqrt(2.))) - tt.sqr(tt.abs_(z)), + tt.log1p(-tt.erfc(z / tt.sqrt(2.))) + ) + class Wald(PositiveContinuous): R""" diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index f06402b4ca..e7c1520544 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -392,12 +392,14 @@ def test_flat(self): def test_normal(self): self.pymc3_matches_scipy(Normal, R, {'mu': R, 'sd': Rplus}, lambda value, mu, sd: sp.norm.logpdf(value, mu, sd)) - self.check_logcdf(Normal, R, {'mu': R, 'sd': Rplus}, + self.check_logcdf(Normal, R, {'mu': R, 'sd': Rplus}, lambda value, mu, sd: sp.norm.logcdf(value, mu, sd)) def test_half_normal(self): self.pymc3_matches_scipy(HalfNormal, Rplus, {'sd': Rplus}, lambda value, sd: sp.halfnorm.logpdf(value, scale=sd)) + self.check_logcdf(HalfNormal, Rplus, {'sd': Rplus}, + lambda value, sd: sp.halfnorm.logcdf(value, scale=sd)) def test_chi_squared(self): self.pymc3_matches_scipy(ChiSquared, Rplus, {'nu': Rplusdunif}, From 14908eb6a03fc99d335afcd2a3202ce6d5f96f21 Mon Sep 17 00:00:00 2001 From: "M. Domenzain" Date: Mon, 6 Mar 2017 20:41:11 +0100 Subject: [PATCH 3/7] Remove useless operation No need to perform both absolute value and square. --- pymc3/distributions/continuous.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index 1d8153d6f7..3d61741fa4 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -252,7 +252,7 @@ def logcdf(self, value): return tt.switch( tt.lt(z, -1.0), tt.log(tt.erfcx(-z / tt.sqrt(2.)) / 2.) - - tt.sqr(tt.abs_(z)) / 2, + tt.sqr(z) / 2, tt.log1p(-tt.erfc(z / tt.sqrt(2.)) / 2.) ) @@ -312,7 +312,7 @@ def logcdf(self, value): z = zvalue(value, mu=0, sd=sd) return tt.switch( tt.lt(z, -1.0), - tt.log(tt.erfcx(-z / tt.sqrt(2.))) - tt.sqr(tt.abs_(z)), + tt.log(tt.erfcx(-z / tt.sqrt(2.))) - tt.sqr(z), tt.log1p(-tt.erfc(z / tt.sqrt(2.))) ) From 46a1cbb9933b1deb420957e73851e6f810e40871 Mon Sep 17 00:00:00 2001 From: "M. Domenzain" Date: Tue, 7 Mar 2017 07:31:43 +0100 Subject: [PATCH 4/7] Implement log CDF for Lognormal distribution --- pymc3/distributions/continuous.py | 16 ++++++++++++++++ pymc3/tests/test_distributions.py | 2 ++ 2 files changed, 18 insertions(+) diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index 3d61741fa4..b068cec687 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -676,6 +676,22 @@ def logp(self, value): - tt.log(value), tau > 0) + def logcdf(self, value): + mu = self.mu + sd = self.sd + z = zvalue(tt.log(value), mu=mu, sd=sd) + + return tt.switch( + tt.le(value, 0), + -np.inf, + tt.switch( + tt.lt(z, -1.0), + tt.log(tt.erfcx(-z / tt.sqrt(2.)) / 2.) - + tt.sqr(z) / 2, + tt.log1p(-tt.erfc(z / tt.sqrt(2.)) / 2.) + ) + ) + class StudentT(Continuous): R""" diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index e7c1520544..01061ec02c 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -461,6 +461,8 @@ def test_lognormal(self): self.pymc3_matches_scipy( Lognormal, Rplus, {'mu': R, 'tau': Rplusbig}, lambda value, mu, tau: sp.lognorm.logpdf(value, tau**-.5, 0, np.exp(mu))) + self.check_logcdf(Lognormal, Rplus, {'mu': R, 'tau': Rplusbig}, + lambda value, mu, tau: sp.lognorm.logcdf(value, tau**-.5, 0, np.exp(mu))) def test_t(self): self.pymc3_matches_scipy(StudentT, R, {'nu': Rplus, 'mu': R, 'lam': Rplus}, From dc8cf206bafe005aab792e2f1dd0ae93cbbe96f2 Mon Sep 17 00:00:00 2001 From: "M. Domenzain" Date: Tue, 7 Mar 2017 07:50:55 +0100 Subject: [PATCH 5/7] Implement log CDF for Cauchy distribution --- pymc3/distributions/continuous.py | 3 +++ pymc3/tests/test_distributions.py | 2 ++ 2 files changed, 5 insertions(+) diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index b068cec687..283c024edc 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -870,6 +870,9 @@ def logp(self, value): - tt.log1p(((value - alpha) / beta)**2), beta > 0) + def logcdf (self, value): + return tt.log(0.5 + tt.arctan ((value - self.alpha) / self.beta) / np.pi) + class HalfCauchy(PositiveContinuous): R""" diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index 01061ec02c..cba4538f09 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -471,6 +471,8 @@ def test_t(self): def test_cauchy(self): self.pymc3_matches_scipy(Cauchy, R, {'alpha': R, 'beta': Rplusbig}, lambda value, alpha, beta: sp.cauchy.logpdf(value, alpha, beta)) + self.check_logcdf(Cauchy, R, {'alpha': R, 'beta': Rplusbig}, + lambda value, alpha, beta: sp.cauchy.logcdf(value, alpha, beta)) def test_half_cauchy(self): self.pymc3_matches_scipy(HalfCauchy, Rplus, {'beta': Rplusbig}, From a9f1b0200515ecf17a1808eb9c1d51095d30f874 Mon Sep 17 00:00:00 2001 From: "M. Domenzain" Date: Tue, 7 Mar 2017 12:40:39 +0100 Subject: [PATCH 6/7] Implement log CDF for Laplace distribution --- pymc3/distributions/continuous.py | 14 ++++++++++++++ pymc3/tests/test_distributions.py | 2 ++ 2 files changed, 16 insertions(+) diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index 283c024edc..4c429e74a3 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -612,6 +612,20 @@ def logp(self, value): return -tt.log(2 * b) - abs(value - mu) / b + def logcdf(self, value): + a = self.mu + b = self.b + y = (value - a) / b + return tt.switch( + tt.le(value, a), + tt.log(0.5) + y, + tt.switch( + tt.gt(y, 1), + tt.log1p(-0.5 * tt.exp(-y)), + tt.log(1 - 0.5 * tt.exp(-y)) + ) + ) + class Lognormal(PositiveContinuous): R""" diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index cba4538f09..54e898a1aa 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -456,6 +456,8 @@ def test_fun(value, mu, alpha): def test_laplace(self): self.pymc3_matches_scipy(Laplace, R, {'mu': R, 'b': Rplus}, lambda value, mu, b: sp.laplace.logpdf(value, mu, b)) + self.check_logcdf(Laplace, R, {'mu': R, 'b': Rplus}, + lambda value, mu, b: sp.laplace.logcdf(value, mu, b)) def test_lognormal(self): self.pymc3_matches_scipy( From 6c76f1fa46dbbdde58744330975b62d2d807705a Mon Sep 17 00:00:00 2001 From: "M. Domenzain" Date: Tue, 7 Mar 2017 14:43:11 +0100 Subject: [PATCH 7/7] Implement log CDF for Triangular distribution --- pymc3/distributions/continuous.py | 18 ++++++++++++++++++ pymc3/tests/test_distributions.py | 2 ++ 2 files changed, 20 insertions(+) diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index 4c429e74a3..206c3ec45d 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -1414,3 +1414,21 @@ def logp(self, value): tt.switch(tt.eq(value, c), tt.log(2 / (upper - lower)), tt.switch(alltrue_elemwise([c < value, value <= upper]), tt.log(2 * (upper - value) / ((upper - lower) * (upper - c))),np.inf))) + + def logcdf(self, value): + l = self.lower + u = self.upper + c = self.c + return tt.switch( + tt.le(value, l), + -np.inf, + tt.switch( + tt.le(value, c), + tt.log(((value - l) ** 2) / ((u - l) * (c - l))), + tt.switch( + tt.lt(value, u), + tt.log1p(-((u - value) ** 2) / ((u - l) * (u - c))), + 0 + ) + ) + ) diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index 54e898a1aa..3e41d821f4 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -373,6 +373,8 @@ def test_triangular(self): self.pymc3_matches_scipy( Triangular, Runif, {'lower': -Rplusunif, 'c': Runif, 'upper': Rplusunif}, lambda value, c, lower, upper: sp.triang.logpdf(value, c-lower, lower, upper-lower)) + self.check_logcdf(Triangular, Runif, {'lower': -Rplusunif, 'c': Runif, 'upper': Rplusunif}, + lambda value, c, lower, upper: sp.triang.logcdf(value, c-lower, lower, upper-lower)) def test_bound_normal(self): PositiveNormal = Bound(Normal, lower=0.)