Scipy.optimize.minimize возвращает неверные результаты

Question

В Python у меня есть функция error_p, которая вычисляет среднеквадратичную ошибку между набором наблюдаемых вероятностей (или, вернее, нормализованными частотами) и распределением Пуассона для данного среднего.

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 

График этой функции (error_p) в диапазоне mu значений:

Очевидно, что есть минимум на входной величины (мю) чуть ниже 2. Однако, когда я называю scipy.optimize.minimize так:

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

Я получаю

(array([ 13.86128699]), 0.007078183160196419) 

, указывающий минимум при мю = 13,86 с функцией величины ~ 0,007, тогда как, если я бегу

error_p(2) 

я получаю

0.000848142902892 

Почему scipy.optimize.minimize не найти истинного минимума?

Answer 1

Если вы используете функцию scipy.optimize.minimize_scalar вы получите ожидаемый результат:

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

Почему scipy.optimize.minimize не работает? Я предполагаю, что ваша функция error_p искажена с точки зрения numpy. Попробуйте следующее:

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

, и вы увидите, что он не работает. Ваша функция не предназначена для ввода массива входов и выплескивает массив выходов, и я думаю, что это то, что ищет минимизация.

Answer 2

Изменить

theory = [poisson.pmf(x, mu) for x in x] 

в

theory = poisson.pmf(x, mu) 

и она работает, как ожидалось.