diff --git a/pymc3/tests/models.py b/pymc3/tests/models.py index a545424d62..6c3d143522 100644 --- a/pymc3/tests/models.py +++ b/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 diff --git a/pymc3/tests/test_examples.py b/pymc3/tests/test_examples.py index 79d3736c79..c0a6fc2ff2 100644 --- a/pymc3/tests/test_examples.py +++ b/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) diff --git a/pymc3/tests/test_starting.py b/pymc3/tests/test_starting.py index 722a16d52c..2e84ba103b 100644 --- a/pymc3/tests/test_starting.py +++ b/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) diff --git a/pymc3/tests/test_tuning.py b/pymc3/tests/test_tuning.py index d4093fbfcd..9a434eab74 100644 --- a/pymc3/tests/test_tuning.py +++ b/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) diff --git a/pymc3/tuning/starting.py b/pymc3/tuning/starting.py index 78a7d64f64..224e7406f5 100644 --- a/pymc3/tuning/starting.py +++ b/pymc3/tuning/starting.py @@ -3,23 +3,27 @@ @author: johnsalvatier ''' -from scipy import optimize 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 -from ..blocking import DictToArrayBijection, ArrayOrdering -from ..util import update_start_vals +from tqdm import tqdm + +from pymc3.vartypes import discrete_types, typefilter +from pymc3.model import modelcontext, Point +from pymc3.blocking import ArrayOrdering, DictToArrayBijection +from pymc3.theanof import inputvars, floatX +from pymc3.util import update_start_vals +from scipy.optimize import minimize + from inspect import getargspec __all__ = ['find_MAP'] -def find_MAP(start=None, vars=None, fmin=None, - return_raw=False, model=None, live_disp=False, callback=None, *args, **kwargs): + +def find_MAP(start=None, vars=None, method=None, progressbar=True, return_raw=False, + model=None, maxeval=50000, callback=None, *args, **kwargs): """ Sets state to the local maximum a posteriori point given a model. Current default of fmin_Hessian does not deal well with optimizing close @@ -30,21 +34,21 @@ 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). - fmin : function - Optimization algorithm (Defaults to `scipy.optimize.fmin_bfgs` unless + method : string or callable + Optimization algorithm (Defaults to `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`) + `Powell` which will perform better). + progressbar : bool + Whether or not to display a progress bar in the command line. + return_raw : bool + Whether to return extra values 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. + maxeval : int + The maximum number of times the posterior distribution is evaluated. 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 + Extra args passed to fmin. """ model = modelcontext(model) if start is None: @@ -69,62 +73,46 @@ def find_MAP(start=None, vars=None, fmin=None, except AttributeError: gradient_avail = False - if disc_vars or not gradient_avail : + if disc_vars or not gradient_avail: pm._log.warning("Warning: gradient not available." + "(E.g. vars contains discrete variables). MAP " + "estimates may not be accurate for the default " + "parameters. Defaulting to non-gradient minimization " + - "fmin_powell.") - fmin = optimize.fmin_powell + "'Powell'.") + method = "Powell" - if fmin is None: + if method is None: if disc_vars: - fmin = optimize.fmin_powell + method = "Powell" else: - fmin = optimize.fmin_bfgs + method = "BFGS" allinmodel(vars, model) 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) - def logp_o(point): - return nan_to_high(-logp(point)) - - # Check to see if minimization function actually uses the gradient - if 'fprime' in getargspec(fmin).args: - dlogp = bij.mapf(model.fastdlogp(vars)) - 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) + if method in ["CG", "BFGS", "Newton-CG", "L-BFGS-B", "TNC", + "SLSQP", "dogleg", "trust-ncg"]: + dlogp_func = bij.mapf(model.fastdlogp(vars)) + cost_func = CostFuncWrapper(maxeval, progressbar, logp_func, dlogp_func) 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) - else: - r = fmin(logp_o, callback=callback, *args, **kwargs) + cost_func = CostFuncWrapper(maxeval, progressbar, logp_func) compute_gradient = False - if isinstance(r, tuple): - mx0 = r[0] - else: - mx0 = r - - if live_disp: - try: - callback.update(mx0) - except: - pass - + try: + r = minimize(cost_func, x0, method=method, jac=compute_gradient, *args, **kwargs) + mx0 = r["x"] + except (KeyboardInterrupt, StopIteration) as e: + mx0, r = cost_func.previous_x, None + cost_func.progress.close() + if isinstance(e, StopIteration): + pm._log.info(e) + finally: + cost_func.progress.close() mx = bij.rmap(mx0) allfinite_mx0 = allfinite(mx0) @@ -169,148 +157,76 @@ 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)) - def nan_to_high(x): return np.where(isfinite(x), x, 1.0e100) - def allinmodel(vars, model): notin = [v for v in vars if v not in model.vars] if notin: 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 = {} +class CostFuncWrapper(object): + def __init__(self, maxeval=5000, progressbar=True, logp_func=None, dlogp_func=None): + self.t0 = time.time() + self.n_eval = 0 + self.maxeval = maxeval + self.logp_func = logp_func + if dlogp_func is None: + self.use_gradient = False + self.desc = 'lp = {:,.5g}' + else: + self.dlogp_func = dlogp_func + self.use_gradient = True + self.desc = 'lp = {:,.5g}, ||grad|| = {:,.5g}' + self.previous_x = None + self.progress = tqdm(total=maxeval, disable=not progressbar) 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""" - - - - - - - - - - - - """.format(self.t_elapsed, self.iters, self.logpost) - if self.dlogp is not None: - html += r""" - - - """.format(self.dlogpost) - html += "
Time Elapsed: {:s}
Iteration: {:d}
Log Posterior: {:.3f}
||grad||: {:.3f}
" - self.prog_table.value = html - ## Parameter table - html = r""" - - - - - - - - - - """ - for var, values in self.paramtable.items(): - html += r""" - - - - - - """.format(var, values["size"], values["valstr"]) - html += "
ParameterSizeCurrent Value
{:s}{:d}{:s}
" - self.param_table.value = html - - -def format_values(val): - fmt = "{:8.3f}" - if val.size == 1: - return fmt.format(np.float(val)) - elif val.size < 9: - return "[" + ", ".join([fmt.format(v) for v in val]) + "]" - else: - start = "[" + ", ".join([fmt.format(v) for v in val[:4]]) - end = ", ".join([fmt.format(v) for v in val[-4:]]) +"]" - return start + ", ... , " + end + neg_value = np.float64(self.logp_func(pm.floatX(x))) + value = -1.0 * nan_to_high(neg_value) + if self.use_gradient: + neg_grad = self.dlogp_func(pm.floatX(x)) + if np.all(np.isfinite(neg_grad)): + self.previous_x = x + grad = nan_to_num(-1.0*neg_grad) + grad = grad.astype(np.float64) + else: + self.previous_x = x + grad = None + + if self.n_eval % 10 == 0: + self.update_progress_desc(neg_value, grad) + + if self.n_eval > self.maxeval: + self.update_progress_desc(neg_value, grad) + self.progress.close() + raise StopIteration + + self.n_eval += 1 + self.progress.update(1) + + if self.use_gradient: + return value, grad + else: + return value + + def update_progress_desc(self, neg_value, grad=None): + if grad is None: + self.progress.set_description(self.desc.format(neg_value)) + else: + norm_grad = np.linalg.norm(grad) + self.progress.set_description(self.desc.format(neg_value, norm_grad)) +