Add MAP estimates for transformed and untransformed variables (#2523)

* Add MAP estimates for transformed and untransformed variables

* Fix docstring

* Fix lint

Author: kyleabeauchamp

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_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,21 @@ 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)
+
+    rtol = 1E-5  # this rtol should work on both floatX precisions
+    np.testing.assert_allclose(map_estimate["mu_i"], truth, rtol=rtol)
+
+    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, rtol=rtol)

pymc3/tuning/starting.py

@@ -7,7 +7,6 @@
 import numpy as np
 from numpy import isfinite, nan_to_num, logical_not
 import pymc3 as pm
-import time
 from ..vartypes import discrete_types, typefilter
 from ..model import modelcontext, Point
 from ..theanof import inputvars
@@ -20,7 +19,7 @@
 
 
 def find_MAP(start=None, vars=None, fmin=None,
-             return_raw=False, model=None, live_disp=False, callback=None,
+             return_raw=False, model=None, callback=None,
              *args, **kwargs):
     """
     Sets state to the local maximum a posteriori point given a model.
@@ -31,20 +30,16 @@ def find_MAP(start=None, vars=None, fmin=None,
     ----------
     start : `dict` of parameter values (Defaults to `model.test_point`)
     vars : list
-        List of variables to set to MAP point (Defaults to all continuous).
+        List of variables to optimize and set to optimum (Defaults to all continuous).
     fmin : function
         Optimization algorithm (Defaults to `scipy.optimize.fmin_bfgs` unless
         discrete variables are specified in `vars`, then
         `scipy.optimize.fmin_powell` which will perform better).
     return_raw : Bool
         Whether to return extra value returned by fmin (Defaults to `False`)
     model : Model (optional if in `with` context)
-    live_disp : Bool
-        Display table tracking optimization progress when run from within
-        an IPython notebook.
     callback : callable
-        Callback function to pass to scipy optimization routine.  Overrides
-        live_disp if callback is given.
+        Callback function to pass to scipy optimization routine.
     *args, **kwargs
         Extra args passed to fmin
     """
@@ -89,7 +84,8 @@ def find_MAP(start=None, vars=None, fmin=None,
 
     start = Point(start, model=model)
     bij = DictToArrayBijection(ArrayOrdering(vars), start)
-
+    logp_func = bij.mapf(model.fastlogp)
+    x0 = bij.map(start)
     logp = bij.mapf(model.fastlogp_nojac)
     def logp_o(point):
         return nan_to_high(-logp(point))
@@ -100,15 +96,9 @@ def logp_o(point):
         def grad_logp_o(point):
             return nan_to_num(-dlogp(point))
 
-        if live_disp and callback is None:
-            callback = Monitor(bij, logp_o, model, grad_logp_o)
-
         r = fmin(logp_o, bij.map(start), fprime=grad_logp_o, callback=callback, *args, **kwargs)
         compute_gradient = True
     else:
-        if live_disp and callback is None:
-            callback = Monitor(bij, logp_o, dlogp=None)
-
         # Check to see if minimization function uses a starting value
         if 'x0' in getargspec(fmin).args:
             r = fmin(logp_o, bij.map(start), callback=callback, *args, **kwargs)
@@ -121,12 +111,6 @@ def grad_logp_o(point):
     else:
         mx0 = r
 
-    if live_disp:
-        try:
-            callback.update(mx0)
-        except:
-            pass
-
     mx = bij.rmap(mx0)
 
     allfinite_mx0 = allfinite(mx0)
@@ -171,13 +155,16 @@ def message(name, values):
                          "density. 2) your distribution logp's are " +
                          "properly specified. Specific issues: \n" +
                          specific_errors)
-    mx = {v.name: mx[v.name].astype(v.dtype) for v in model.vars}
+
+    vars = model.unobserved_RVs
+    mx = {var.name: value for var, value in zip(vars, model.fastfn(vars)(mx))}
 
     if return_raw:
         return mx, r
     else:
         return mx
 
+
 def allfinite(x):
     return np.all(isfinite(x))
 
@@ -192,120 +179,6 @@ def allinmodel(vars, model):
         raise ValueError("Some variables not in the model: " + str(notin))
 
 
-
-class Monitor(object):
-    def __init__(self, bij, logp, model, dlogp=None):
-        try:
-            from IPython.display import display
-            from ipywidgets import HTML, VBox, HBox, FlexBox
-            self.prog_table  = HTML(width='100%')
-            self.param_table = HTML(width='100%')
-            r_col = VBox(children=[self.param_table], padding=3, width='100%')
-            l_col = HBox(children=[self.prog_table],  padding=3, width='25%')
-            self.hor_align = FlexBox(children = [l_col, r_col], width='100%', orientation='vertical')
-            display(self.hor_align)
-            self.using_notebook = True
-            self.update_interval = 1
-        except:
-            self.using_notebook = False
-            self.update_interval = 2
-
-        self.iters = 0
-        self.bij = bij
-        self.model = model
-        self.fn = model.fastfn(model.unobserved_RVs)
-        self.logp = logp
-        self.dlogp = dlogp
-        self.t_initial = time.time()
-        self.t0 = self.t_initial
-        self.paramtable = {}
-
-    def __call__(self, x):
-        self.iters += 1
-        if time.time() - self.t0 > self.update_interval or self.iters == 1:
-            self.update(x)
-
-    def update(self, x):
-        self._update_progtable(x)
-        self._update_paramtable(x)
-        if self.using_notebook:
-            self._display_notebook()
-            self.t0 = time.time()
-
-    def _update_progtable(self, x):
-        s = time.time() - self.t_initial
-        hours, remainder = divmod(int(s), 3600)
-        minutes, seconds = divmod(remainder, 60)
-        self.t_elapsed = "{:2d}h{:2d}m{:2d}s".format(hours, minutes, seconds)
-        self.logpost = -1.0*np.float(self.logp(x))
-        self.dlogpost = np.linalg.norm(self.dlogp(x))
-
-    def _update_paramtable(self, x):
-        var_state = self.fn(self.bij.rmap(x))
-        for var, val in zip(self.model.unobserved_RVs, var_state):
-            if not var.name.endswith("_"):
-                valstr = format_values(val)
-                self.paramtable[var.name] = {"size": val.size, "valstr": valstr}
-
-    def _display_notebook(self):
-        ## Progress table
-        html = r"""<style type="text/css">
-        table { border-collapse:collapse }
-        .tg {border-collapse:collapse;border-spacing:0;border:none;}
-        .tg td{font-family:Arial, sans-serif;font-size:14px;padding:3px 3px;border-style:solid;border-width:0px;overflow:hidden;word-break:normal;}
-        .tg th{Impact, Charcoal, sans-serif;font-size:13px;font-weight:bold;padding:3px 3px;border-style:solid;border-width:0px;overflow:hidden;word-break:normal; background-color:#0E688A;color:#ffffff;}
-        .tg .tg-vkoh{white-space:pre;font-weight:normal;font-family:"Lucida Console", Monaco, monospace !important; background-color:#ffffff;color:#000000}
-        .tg .tg-suao{font-weight:bold;font-family:"Lucida Console", Monaco, monospace !important;background-color:#0E688A;color:#ffffff;}
-        """
-        html += r"""
-        </style>
-        <table class="tg" style="undefined;">
-           <col width="400px" />
-           <tr>
-             <th class= "tg-vkoh">Time Elapsed: {:s}</th>
-           </tr>
-           <tr>
-             <th class= "tg-vkoh">Iteration: {:d}</th>
-           </tr>
-           <tr>
-             <th class= "tg-vkoh">Log Posterior: {:.3f}</th>
-           </tr>
-        """.format(self.t_elapsed, self.iters, self.logpost)
-        if self.dlogp is not None:
-           html += r"""
-             <tr>
-               <th class= "tg-vkoh">||grad||: {:.3f}</th>
-             </tr>""".format(self.dlogpost)
-        html += "</table>"
-        self.prog_table.value = html
-        ## Parameter table
-        html = r"""<style type="text/css">
-          .tg .tg-bgft{font-weight:normal;font-family:"Lucida Console", Monaco, monospace !important;background-color:#0E688A;color:#ffffff;}
-          .tg td{font-family:Arial, sans-serif;font-size:12px;padding:3px 3px;border-style:solid;border-width:1px;overflow:hidden;word-break:normal;border-color:#504A4E;color:#333;background-color:#fff;word-wrap: break-word;}
-          .tg th{Impact, Charcoal, sans-serif;font-size:13px;font-weight:bold;padding:3px 3px;border-style:solid;border-width:1px;overflow:hidden;word-break:normal;border-color:#504A4E;background-color:#0E688A;color:#ffffff;}
-          </style>
-          <table class="tg" style="undefined;">
-             <col width="130px" />
-             <col width="50px" />
-             <col width="600px" />
-             <tr>
-               <th class="tg">Parameter</th>
-               <th class="tg">Size</th>
-               <th class="tg">Current Value</th>
-             </tr>
-           """
-        for var, values in self.paramtable.items():
-            html += r"""
-              <tr>
-                <td class="tg-bgft">{:s}</td>
-                <td class="tg-vkoh">{:d}</td>
-                <td class="tg-vkoh">{:s}</td>
-              </tr>
-            """.format(var, values["size"], values["valstr"])
-        html += "</table>"
-        self.param_table.value = html
-
-
 def format_values(val):
     fmt = "{:8.3f}"
     if val.size == 1: