Merge branch 'master' into sample_arg_rename

Author: fonnesbeck

RELEASE-NOTES.md

@@ -8,6 +8,11 @@
   of the success probability. This is faster and more stable than using
   `p=tt.nnet.sigmoid(logit_p)`.
 
+### Fixes
+
+- `VonMises` does not overflow for large values of kappa. i0 and i1 have been removed and we now use
+   log_i0 to compute the logp.
+
 ### Deprecations
 
 - DIC and BPIC calculations have been removed

pymc3/backends/report.py

@@ -1,6 +1,7 @@
 from collections import namedtuple
 import logging
 import enum
+from ..util import is_transformed_name, get_untransformed_name
 
 
 logger = logging.getLogger('pymc3')
@@ -69,10 +70,15 @@ def _run_convergence_checks(self, trace, model):
 
         from pymc3 import diagnostics
 
-        varname = [rv.name for rv in model.free_RVs]
+        varnames = []
+        for rv in model.free_RVs:
+            rv_name = rv.name
+            if is_transformed_name(rv_name):
+                rv_name = get_untransformed_name(rv_name)
+            varnames.append(rv_name)
 
-        self._effective_n = effective_n = diagnostics.effective_n(trace, varname)
-        self._gelman_rubin = gelman_rubin = diagnostics.gelman_rubin(trace, varname)
+        self._effective_n = effective_n = diagnostics.effective_n(trace, varnames)
+        self._gelman_rubin = gelman_rubin = diagnostics.gelman_rubin(trace, varnames)
 
         warnings = []
         rhat_max = max(val.max() for val in gelman_rubin.values())

pymc3/distributions/continuous.py

@@ -16,7 +16,7 @@
 from pymc3.theanof import floatX
 from . import transforms
 from pymc3.util import get_variable_name
-from .special import i0, i1
+from .special import log_i0
 from .dist_math import bound, logpow, gammaln, betaln, std_cdf, alltrue_elemwise, SplineWrapper
 from .distribution import Continuous, draw_values, generate_samples
 
@@ -1850,8 +1850,7 @@ def __init__(self, mu=0.0, kappa=None, transform='circular',
             transform = transforms.Circular()
         super(VonMises, self).__init__(transform=transform, *args, **kwargs)
         self.mean = self.median = self.mode = self.mu = mu = tt.as_tensor_variable(mu)
-        self.kappa = kappa = tt.as_tensor_variable(kappa)
-        self.variance = 1 - i1(kappa) / i0(kappa)
+        self.kappa = kappa = floatX(tt.as_tensor_variable(kappa))
 
         assert_negative_support(kappa, 'kappa', 'VonMises')
 
@@ -1865,7 +1864,8 @@ def random(self, point=None, size=None, repeat=None):
     def logp(self, value):
         mu = self.mu
         kappa = self.kappa
-        return bound(kappa * tt.cos(mu - value) - tt.log(2 * np.pi * i0(kappa)), value >= -np.pi, value <= np.pi, kappa >= 0)
+        return bound(kappa * tt.cos(mu - value) - (tt.log(2 * np.pi) + log_i0(kappa)),
+                     kappa > 0, value >= -np.pi, value <= np.pi)
 
     def _repr_latex_(self, name=None, dist=None):
         if dist is None:

pymc3/distributions/special.py

@@ -1,9 +1,9 @@
 import numpy as np
 import theano.tensor as tt
-from theano.scalar.basic_scipy import GammaLn, Psi, I0, I1
+from theano.scalar.basic_scipy import GammaLn, Psi
 from theano import scalar
 
-__all__ = ['gammaln', 'multigammaln', 'psi', 'i0', 'i1']
+__all__ = ['gammaln', 'multigammaln', 'psi', 'log_i0']
 
 scalar_gammaln = GammaLn(scalar.upgrade_to_float, name='scalar_gammaln')
 gammaln = tt.Elemwise(scalar_gammaln, name='gammaln')
@@ -12,20 +12,28 @@
 def multigammaln(a, p):
     """Multivariate Log Gamma
 
-    :Parameters:
-        a : tensor like
-        p : int degrees of freedom
-            p > 0
+    Parameters
+    ----------
+    a : tensor like
+    p : int
+       degrees of freedom. p > 0
     """
     i = tt.arange(1, p + 1)
     return (p * (p - 1) * tt.log(np.pi) / 4.
             + tt.sum(gammaln(a + (1. - i) / 2.), axis=0))
 
-scalar_psi = Psi(scalar.upgrade_to_float, name='scalar_psi')
-psi = tt.Elemwise(scalar_psi, name='psi')
 
-scalar_i0 = I0(scalar.upgrade_to_float, name='scalar_i0')
-i0 = tt.Elemwise(scalar_i0, name='i0')
+def log_i0(x):
+    """
+    Calculates the logarithm of the 0 order modified Bessel function of the first kind""
+    """
+    return tt.switch(tt.lt(x, 5), tt.log1p(x**2. / 4. + x**4. / 64. + x**6. / 2304.
+                                           + x**8. / 147456. + x**10. / 14745600.
+                                           + x**12. / 2123366400.),
+                                  x - 0.5 * tt.log(2. * np.pi * x) + tt.log1p(1. / (8. * x)
+                                  + 9. / (128. * x**2.) + 225. / (3072. * x**3.)
+                                  + 11025. / (98304. * x**4.)))
 
-scalar_i1 = I1(scalar.upgrade_to_float, name='scalar_i1')
-i1 = tt.Elemwise(scalar_i1, name='i1')
+
+scalar_psi = Psi(scalar.upgrade_to_float, name='scalar_psi')
+psi = tt.Elemwise(scalar_psi, name='psi')

pymc3/sampling.py

@@ -258,8 +258,9 @@ def sample(draws=500, step=None, init='auto', n_init=200000, start=None,
     tune : int
         Number of iterations to tune, if applicable (defaults to 500).
         Samplers adjust the step sizes, scalings or similar during
-        tuning. These samples will be drawn in addition to samples
-        and discarded unless discard_tuned_samples is set to True.
+        tuning. Tuning samples will be drawn in addition to the number 
+        specified in the `draws` argument, and will be discarded 
+        unless `discard_tuned_samples` is set to False.
     nuts_kwargs : dict
         Options for the NUTS sampler. See the docstring of NUTS
         for a complete list of options. Common options are

pymc3/tests/test_distributions.py

@@ -931,6 +931,7 @@ def test_ex_gaussian(self, value, mu, sigma, nu, logp):
         pt = {'eg': value}
         assert_almost_equal(model.fastlogp(pt), logp, decimal=select_by_precision(float64=6, float32=2), err_msg=str(pt))
 
+    @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32")
     def test_vonmises(self):
         self.pymc3_matches_scipy(
             VonMises, R, {'mu': Circ, 'kappa': Rplus},