Scipy Error when computing MAP

[plison]

I get an error when computing MAP estimates:

import pymc3 as pm
with pm.Model() as model:
    rain = pm.Bernoulli('rain', 0.2)
    pm.find_MAP()

which throws the following error in Scipy:

error                                     Traceback (most recent call last)
<ipython-input-1-d6fe66350f1d> in <module>()
      3 with pm.Model() as model:
      4     rain = pm.Bernoulli('rain', 0.2)
----> 5     pm.find_MAP()

~/anaconda3/lib/python3.6/site-packages/pymc3/tuning/starting.py in find_MAP(start, vars, method, return_raw, include_transformed, progressbar, maxeval, model, *args, **kwargs)
    120 
    121         try:
--> 122             opt_result = minimize(cost_func, x0, method=method, jac=compute_gradient, *args, **kwargs)
    123             mx0 = opt_result["x"]  # r -> opt_result
    124             cost_func.progress.total = cost_func.progress.n + 1

~/anaconda3/lib/python3.6/site-packages/scipy/optimize/_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    601     elif meth == 'l-bfgs-b':
    602         return _minimize_lbfgsb(fun, x0, args, jac, bounds,
--> 603                                 callback=callback, **options)
    604     elif meth == 'tnc':
    605         return _minimize_tnc(fun, x0, args, jac, bounds, callback=callback,

~/anaconda3/lib/python3.6/site-packages/scipy/optimize/lbfgsb.py in _minimize_lbfgsb(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, **unknown_options)
    326         _lbfgsb.setulb(m, x, low_bnd, upper_bnd, nbd, f, g, factr,
    327                        pgtol, wa, iwa, task, iprint, csave, lsave,
--> 328                        isave, dsave, maxls)
    329         task_str = task.tostring()
    330         if task_str.startswith(b'FG'):

error: failed in converting 3rd argument `l' of _lbfgsb.setulb to C/Fortran array

I'm currently running the latest development version of PyMC3 (but I also tried the stable version on anaconda, and I got the same error). The scipy version is 1.1.0.

Any clue as to why this is happening?

[aseyboldt]

This is because you are using find_MAP on a model with unobserved discrete variables. We use the scipy optimisers internally, and they assume that the optimisation problem is continuous, so we don't support that (discrete parameters are also typically problematic for sampling by the way. Often there are ways around using them.) We should really give a better error message here though...

[aseyboldt]

@bwengals I think you worked on this last, right? It seems there is some code here that should switch to method="Powell" if there are discrete parameters (this code somehow doesn't work). This still isn't the right thing though. Powell doesn't use the gradient, but it still assumes that the variables are continuous:

optimize.minimize(fun=lambda x: (x - 0.1)**2, x0=np.array([3]), method='Powell')
   direc: array([[1.]])
     fun: array(2.86655413e-30)
 message: 'Optimization terminated successfully.'
    nfev: 20
     nit: 2
  status: 0
 success: True
       x: array(0.1)

Should we just throw an error instead or at least print a warning? Can continuous values still be useful in some settings?