Change tune method to not act inplace plus more sampler stats (#3732)

* add regression test to demo #3731

* return new array and don't touch original in tune method
+ see #3731

* add 'accepted' and 'scaling'  to (DE)Metropolis sampler stats

* use existing model.ndim property for code clarity

Author: michaelosthege

pymc3/step_methods/metropolis.py

@@ -93,7 +93,9 @@ class Metropolis(ArrayStepShared):
     generates_stats = True
     stats_dtypes = [{
         'accept': np.float64,
+        'accepted': np.bool,
         'tune': np.bool,
+        'scaling': np.float64,
     }]
 
     def __init__(self, vars=None, S=None, proposal_dist=None, scaling=1.,
@@ -166,7 +168,9 @@ def astep(self, q0):
 
         stats = {
             'tune': self.tune,
+            'scaling': self.scaling,
             'accept': np.exp(accept),
+            'accepted': accepted,
         }
 
         return q_new, [stats]
@@ -191,26 +195,24 @@ def tune(scale, acc_rate):
     >0.95         x 10
 
     """
-
-    # Switch statement
     if acc_rate < 0.001:
         # reduce by 90 percent
-        scale *= 0.1
+        return scale * 0.1
     elif acc_rate < 0.05:
         # reduce by 50 percent
-        scale *= 0.5
+        return scale * 0.5
     elif acc_rate < 0.2:
         # reduce by ten percent
-        scale *= 0.9
+        return scale * 0.9
     elif acc_rate > 0.95:
         # increase by factor of ten
-        scale *= 10.0
+        return scale * 10.0
     elif acc_rate > 0.75:
         # increase by double
-        scale *= 2.0
+        return scale * 2.0
     elif acc_rate > 0.5:
         # increase by ten percent
-        scale *= 1.1
+        return scale * 1.1
 
     return scale
 
@@ -531,7 +533,9 @@ class DEMetropolis(PopulationArrayStepShared):
     generates_stats = True
     stats_dtypes = [{
         'accept': np.float64,
+        'accepted': np.bool,
         'tune': np.bool,
+        'scaling': np.float64,
     }]
 
     def __init__(self, vars=None, S=None, proposal_dist=None, lamb=None, scaling=0.001,
@@ -544,7 +548,7 @@ def __init__(self, vars=None, S=None, proposal_dist=None, lamb=None, scaling=0.0
         vars = pm.inputvars(vars)
 
         if S is None:
-            S = np.ones(sum(v.dsize for v in vars))
+            S = np.ones(model.ndim)
 
         if proposal_dist is not None:
             self.proposal_dist = proposal_dist(S)
@@ -553,7 +557,7 @@ def __init__(self, vars=None, S=None, proposal_dist=None, lamb=None, scaling=0.0
 
         self.scaling = np.atleast_1d(scaling).astype('d')
         if lamb is None:
-            lamb = 2.38 / np.sqrt(2 * S.size)
+            lamb = 2.38 / np.sqrt(2 * model.ndim)
         self.lamb = float(lamb)
         self.tune = tune
         self.tune_interval = tune_interval
@@ -593,7 +597,9 @@ def astep(self, q0):
 
         stats = {
             'tune': self.tune,
+            'scaling': self.scaling,
             'accept': np.exp(accept),
+            'accepted': accepted
         }
 
         return q_new, [stats]

pymc3/tests/test_tuning.py

@@ -1,5 +1,6 @@
 import numpy as np
 from numpy import inf
+from pymc3.step_methods.metropolis import tune
 from pymc3.tuning import scaling, find_MAP
 from . import models
 
@@ -32,3 +33,11 @@ def test_mle_jacobian():
         map_estimate = find_MAP(method="BFGS", model=model)
 
     np.testing.assert_allclose(map_estimate["mu_i"], truth, rtol=rtol)
+
+
+def test_tune_not_inplace():
+    orig_scaling = np.array([0.001, 0.1])
+    returned_scaling = tune(orig_scaling, acc_rate=0.6)
+    assert not returned_scaling is orig_scaling
+    assert np.all(orig_scaling == np.array([0.001, 0.1]))
+    pass