From f942f1e18bfe8c75854cca59148294c642ea5467 Mon Sep 17 00:00:00 2001 From: Chris Fonnesbeck Date: Tue, 18 Oct 2016 14:51:54 -0500 Subject: [PATCH 1/4] More robust transformations and scaling adjustments --- pymc3/distributions/continuous.py | 3 +-- pymc3/distributions/transforms.py | 9 +++++---- pymc3/tuning/scaling.py | 14 +++++++------- 3 files changed, 13 insertions(+), 13 deletions(-) diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index bf6ce60ab5..0294fb2782 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -217,10 +217,9 @@ def random(self, point=None, size=None, repeat=None): def logp(self, value): tau = self.tau - sd = self.sd mu = self.mu return bound((-tau * (value - mu)**2 + tt.log(tau / np.pi / 2.)) / 2., - tau > 0, sd > 0) + tau > 0) class HalfNormal(PositiveContinuous): diff --git a/pymc3/distributions/transforms.py b/pymc3/distributions/transforms.py index 00b95b6582..ad2a063d5c 100644 --- a/pymc3/distributions/transforms.py +++ b/pymc3/distributions/transforms.py @@ -109,18 +109,19 @@ class Interval(ElemwiseTransform): name = "interval" - def __init__(self, a, b): + def __init__(self, a, b, eps=1e-6): self.a = a self.b = b + self.eps = eps def backward(self, x): a, b = self.a, self.b - r = (b - a) * tt.exp(x) / (1 + tt.exp(x)) + a + r = (b - a) / (1 + tt.exp(-x)) + a return r def forward(self, x): - a, b = self.a, self.b - r = tt.log((x - a) / (b - x)) + a, b, e = self.a, self.b, self.eps + r = tt.log(tt.maximum((x - a) / tt.maximum(b - x, e), e)) return r interval = Interval diff --git a/pymc3/tuning/scaling.py b/pymc3/tuning/scaling.py index 654b8f5aca..5e37c0412b 100644 --- a/pymc3/tuning/scaling.py +++ b/pymc3/tuning/scaling.py @@ -105,28 +105,28 @@ def find_hessian_diag(point, vars=None, model=None): return H(Point(point, model=model)) -def guess_scaling(point, vars=None, model=None): +def guess_scaling(point, vars=None, model=None, scaling_bound=1e-3): model = modelcontext(model) try: h = find_hessian_diag(point, vars, model=model) except NotImplementedError: h = fixed_hessian(point, vars, model=model) - return adjust_scaling(h) + return adjust_scaling(h, scaling_bound) -def adjust_scaling(s): +def adjust_scaling(s, bound): if s.ndim < 2: - return adjust_precision(s) + return adjust_precision(s, bound) else: val, vec = np.linalg.eigh(s) - val = adjust_precision(val) + val = adjust_precision(val, bound) return eig_recompose(val, vec) -def adjust_precision(tau): +def adjust_precision(tau, bound): mag = sqrt(abs(tau)) - bounded = bound(log(mag), log(1e-10), log(1e10)) + bounded = bound(log(mag), log(1/bound), log(bound)) return exp(bounded)**2 From c8dc21f6bc5fe0b8359817ce063caa073ce612e5 Mon Sep 17 00:00:00 2001 From: Chris Fonnesbeck Date: Tue, 18 Oct 2016 15:19:21 -0500 Subject: [PATCH 2/4] Fixed bad name for scaling bound --- pymc3/tuning/scaling.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pymc3/tuning/scaling.py b/pymc3/tuning/scaling.py index 5e37c0412b..3a92375b04 100644 --- a/pymc3/tuning/scaling.py +++ b/pymc3/tuning/scaling.py @@ -114,19 +114,19 @@ def guess_scaling(point, vars=None, model=None, scaling_bound=1e-3): return adjust_scaling(h, scaling_bound) -def adjust_scaling(s, bound): +def adjust_scaling(s, scaling_bound): if s.ndim < 2: - return adjust_precision(s, bound) + return adjust_precision(s, scaling_bound) else: val, vec = np.linalg.eigh(s) - val = adjust_precision(val, bound) + val = adjust_precision(val, scaling_bound) return eig_recompose(val, vec) -def adjust_precision(tau, bound): +def adjust_precision(tau, scaling_bound): mag = sqrt(abs(tau)) - bounded = bound(log(mag), log(1/bound), log(bound)) + bounded = bound(log(mag), log(1/scaling_bound), log(scaling_bound)) return exp(bounded)**2 From 1a650aa2a24639ea16c453672bd99a14b32901b8 Mon Sep 17 00:00:00 2001 From: Chris Fonnesbeck Date: Tue, 18 Oct 2016 15:53:33 -0500 Subject: [PATCH 3/4] Tests sd rather than tau (or both) in normal logp --- pymc3/distributions/continuous.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index 0294fb2782..11df8121ec 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -216,10 +216,11 @@ def random(self, point=None, size=None, repeat=None): size=size) def logp(self, value): + sd = self.sd tau = self.tau mu = self.mu return bound((-tau * (value - mu)**2 + tt.log(tau / np.pi / 2.)) / 2., - tau > 0) + sd > 0) class HalfNormal(PositiveContinuous): From 4450043b7ec107731f071c699c1ab56206c99690 Mon Sep 17 00:00:00 2001 From: Chris Fonnesbeck Date: Tue, 18 Oct 2016 17:52:45 -0500 Subject: [PATCH 4/4] Added default scaling_bound for adjust_precision --- pymc3/tuning/scaling.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pymc3/tuning/scaling.py b/pymc3/tuning/scaling.py index 3a92375b04..736ccfd316 100644 --- a/pymc3/tuning/scaling.py +++ b/pymc3/tuning/scaling.py @@ -123,7 +123,7 @@ def adjust_scaling(s, scaling_bound): return eig_recompose(val, vec) -def adjust_precision(tau, scaling_bound): +def adjust_precision(tau, scaling_bound=1e-3): mag = sqrt(abs(tau)) bounded = bound(log(mag), log(1/scaling_bound), log(scaling_bound))