Scipy.optimize.minimize returning incorrect results

Question

In Python, I have a function error_p that calculates the mean squared error between a set of observed probabilities (or more correctly, normalised frequencies), and the Poisson distribution for a given mean.

from scipy.stats import poisson
import numpy as np
from scipy.optimize import minimize

def error_p(mu):
    """ 
    Returns mean squared error between observed data points and Poisson distribution 
    """
    data = np.array([17.0,32,20,19,6,5,7,5,0,1,3,1,1,0,0,0,0,1,0,0])
    data = data / sum(data)
    x = range(len(data))
    theory = [poisson.pmf(x, mu) for x in x]
    error = np.mean((data - theory)**2)
    return error

A plot of this function (error_p) over a range of mu values is:

Clearly there's a minimum at an input (mu) value of just below 2. However, when I call scipy.optimize.minimize like this:

results = minimize(error_p, 2, tol=0.00001)
results['x'], results['fun']

I get

(array([ 13.86128699]), 0.007078183160196419)

indicating a minimum at mu=13.86 with a function value of ~ 0.007, whereas if I run

error_p(2)

I get

0.000848142902892

Why is scipy.optimize.minimize not finding the true minimum?

Answer 1

If you use the function scipy.optimize.minimize_scalar you get the expected result:

results = minimize_scalar(error_p, tol=0.00001)
print results['x'], results['fun']
>>> 1.88536329298 0.000820148069544

Why does scipy.optimize.minimize not work? My guess is that your function error_p is malformed from a numpy perspective. Try this:

MU = np.linspace(0,20,100)
error_p(MU)

and you'll see that it fails. Your function isn't tailored to take in an array of inputs and spit out an array of outputs, which I think is what minimize is looking for.