Максимизация функции

Question

Я пытаюсь максимизировать функцию, используя следующий код:

from scipy.optimize import minimize
from scipy.stats import lognorm, norm
import numpy as np

np.random.seed(123)
obs = np.random.normal(loc=20, scale=3, size=20)

# Log-Posterior optimisation objectiv 
def objective(params, y):
    mu = params[0]
    sigma = params[1]
    llikelihood = np.sum(np.log(norm.pdf(y, mu, sigma)))
    lpost = llikelihood + np.log(norm.pdf(mu, 0, 100)) + np.log(lognorm.pdf(sigma, loc= 0, s = 4))
    return -1*lpost

starting_mu = 0
starting_sigma = 1
optim_res = minimize(fun = objective, x0=(starting_mu, starting_sigma), args=(obs))

Код работает до последней строки оптимизации.Я совершенно уверен, что ошибка заключается в том, что я пытаюсь выполнить оптимизацию, как в R, используя идентичные настройки и наблюдения, objective() оценивает одно и то же значение.Кроме того, при использовании optim() функция оптимизируется до значения mu = 21,6 и sigma = 3,28.

Я могу использовать код R, однако было бы проще запустить код на Python, чтобы он мог интегрироваться со всем остальным, что я делаю.

РЕДАКТИРОВАТЬ: сообщение трассировки:

dert2@ma0phd201803:~$ python laplace_approx.py 
laplace_approx.py:12: RuntimeWarning: divide by zero encountered in log
    lpost = llikelihood + np.log(norm.pdf(mu, 0, 100)) + np.log(lognorm.pdf(sigma, loc= 0, s = 4))
laplace_approx.py:12: RuntimeWarning: divide by zero encountered in log
    lpost = llikelihood + np.log(norm.pdf(mu, 0, 100)) + np.log(lognorm.pdf(sigma, loc= 0, s = 4))
laplace_approx.py:12: RuntimeWarning: divide by zero encountered in log
    lpost = llikelihood + np.log(norm.pdf(mu, 0, 100)) + np.log(lognorm.pdf(sigma, loc= 0, s = 4))
/opt/anaconda3/lib/python3.7/site-packages/numpy/core/fromnumeric.py:83:     RuntimeWarning: invalid value encountered in reduce
    return ufunc.reduce(obj, axis, dtype, out, **passkwargs)

Answer 1

Некоторые быстрые шаги в отладчике показывают, что norm.pdf(y, mu, sigma) возвращает массив nan за одну итерацию.На этой же итерации lognorm.pdf(sigma, loc=0, s=4) возвращает ноль.Еще до этого sigma было присвоено отрицательное значение.Это действительно, где проблема начинается.Потребовалось бы некоторое усилие, чтобы точно отследить, почему функция вызывается с отрицательным значением params[1], но это происходит.Я бы использовал отладчик, чтобы пройти через функцию minimize, чтобы увидеть, что происходит.Я считаю, что отрицательное значение передается при 9-м вызове функции, если это помогает.Возможно, что метод, который используется для оптимизации функции, не подходит для вашей конкретной функции.Вы можете попробовать использовать другой метод, выбранный путем присвоения параметра method (подробности см. В документации).