diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py
index bf6ce60ab5..11df8121ec 100644
--- a/pymc3/distributions/continuous.py
+++ b/pymc3/distributions/continuous.py
@@ -216,11 +216,11 @@ def random(self, point=None, size=None, repeat=None):
                                 size=size)
 
     def logp(self, value):
-        tau = self.tau
         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)
+                     sd > 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..736ccfd316 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, scaling_bound):
     if s.ndim < 2:
-        return adjust_precision(s)
+        return adjust_precision(s, scaling_bound)
     else:
         val, vec = np.linalg.eigh(s)
-        val = adjust_precision(val)
+        val = adjust_precision(val, scaling_bound)
         return eig_recompose(val, vec)
 
 
-def adjust_precision(tau):
+def adjust_precision(tau, scaling_bound=1e-3):
     mag = sqrt(abs(tau))
 
-    bounded = bound(log(mag), log(1e-10), log(1e10))
+    bounded = bound(log(mag), log(1/scaling_bound), log(scaling_bound))
     return exp(bounded)**2