Add test for MAP/MLE estimation

pymc3/tests/models.py

@@ -157,3 +157,19 @@ def beta_bernoulli(n=2):
         pm.Beta('x', 3, 1, shape=n, transform=None)
         pm.Bernoulli('y', 0.5)
     return model.test_point, model, None
+
+
+def simple_normal(bounded_prior=False):
+    """Simple normal for testing MLE / MAP; probes issue #2482."""
+    x0 = 10.0
+    sd = 1.0
+    a, b = (9, 12)  # bounds for uniform RV, need non-symmetric to reproduce issue
+
+    with pm.Model() as model:
+        if bounded_prior:
+            mu_i = pm.Uniform("mu_i", a, b)
+        else:
+            mu_i = pm.Flat("mu_i")
+        pm.Normal("X_obs", mu=mu_i, sd=sd, observed=x0)
+
+    return model.test_point, model, None

pymc3/tests/test_examples.py

@@ -195,7 +195,7 @@ def build_model(self):
 
     def test_run(self):
         with self.build_model():
-            start = pm.find_MAP(fmin=opt.fmin_powell)
+            start = pm.find_MAP(method="powell")
             pm.sample(50, pm.Slice(), start=start)
 
 

pymc3/tests/test_starting.py

@@ -75,7 +75,7 @@ def test_find_MAP():
         # Test gradient minimization
         map_est1 = starting.find_MAP()
         # Test non-gradient minimization
-        map_est2 = starting.find_MAP(fmin=starting.optimize.fmin_powell)
+        map_est2 = starting.find_MAP(method="powell")
 
     close_to(map_est1['mu'], 0, tol)
     close_to(map_est1['sigma'], 1, tol)

pymc3/tests/test_tuning.py

@@ -1,6 +1,6 @@
 import numpy as np
 from numpy import inf
-from pymc3.tuning import scaling
+from pymc3.tuning import scaling, find_MAP
 from . import models
 
 
@@ -14,3 +14,20 @@ def test_guess_scaling():
     start, model, _ = models.non_normal(n=5)
     a1 = scaling.guess_scaling(start, model=model)
     assert all((a1 > 0) & (a1 < 1e200))
+
+
+def test_mle_jacobian():
+    """Test MAP / MLE estimation for distributions with flat priors."""
+    truth = 10.0  # Simple normal model should give mu=10.0
+
+    start, model, _ = models.simple_normal(bounded_prior=False)
+    with model:
+        map_estimate = find_MAP(model=model)
+
+    np.testing.assert_allclose(map_estimate["mu_i"], truth)
+
+    start, model, _ = models.simple_normal(bounded_prior=True)
+    with model:
+        map_estimate = find_MAP(model=model)
+
+    np.testing.assert_allclose(map_estimate["mu_i"], truth)