diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index 746d4aa240..206c3ec45d 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): """ @@ -240,10 +252,10 @@ 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.) ) - + class HalfNormal(PositiveContinuous): R""" @@ -295,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(z), + tt.log1p(-tt.erfc(z / tt.sqrt(2.))) + ) + class Wald(PositiveContinuous): R""" @@ -591,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""" @@ -655,6 +690,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""" @@ -833,6 +884,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""" @@ -1360,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 9449bb8a78..3e41d821f4 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -366,11 +366,15 @@ 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( 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.) @@ -390,12 +394,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}, @@ -452,11 +458,15 @@ 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( 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}, @@ -465,6 +475,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},