Как обеспечить числовую стабильность при минимизации (NLL)? - PullRequest
Новая
       6

Как обеспечить числовую стабильность при минимизации (NLL)?

0 голосов
/ 26 августа 2018

Если мы возьмем пример scipy.optimize.curve_fit и слегка изменим его, чтобы параметры наилучшего соответствия были выбраны в качестве Оценщиков максимального правдоподобия (MLE) и были найдены с использованием scipy.optimize.minimizeс функцией потерь, представляющей отрицательное логарифмическое правдоподобие, кажется, есть проблемы с числовой стабильностью: 

import numpy as np
from scipy.optimize import curve_fit
from scipy.optimize import minimize

np.random.seed(1729)


def func(x, a, b, c):
    return a * np.exp(-b * x) + c


def NLL(p, data):
    log_likelihood = np.vectorize(
        lambda x: np.log(p[0] * np.exp(-p[1] * x) + p[2]))
    return -1. * np.array(log_likelihood(data)).sum()


def main():
    x_data = np.linspace(0, 4, 50)
    y = func(x_data, 2.5, 1.3, 0.5)
    y_noise = 0.2 * np.random.normal(size=x_data.size)
    y_data = y + y_noise

    init_params = [3., 1.5, 0.5]
    print('\n### Using minimize\n')
    minimize_result = minimize(NLL, x0=init_params, args=(x_data),
                               method='BFGS', options={'disp': True})
    print('\n')
    print(minimize_result)
    print('\n### Using curve_fit\n')

    init_params = [3., 1.5, 0.5]
    popt, pcov = curve_fit(func, x_data, y_data,
                           p0=init_params, bounds=(0, [4., 2., 0.5]))
    print('fit values: {}'.format(popt))
    print('covariance matrix:\n{}'.format(pcov))
    print('uncertainties: {}\n'.format(np.sqrt(np.diag(pcov))))


if __name__ == '__main__':
    main()

производит 

### Using minimize

/home/mcf/anaconda3/lib/python3.5/site-packages/scipy/optimize/optimize.py:663: RuntimeWarning: invalid value encountered in double_scalars
  grad[k] = (f(*((xk + d,) + args)) - f0) / d[k]
curve_fitting_example.py:16: RuntimeWarning: overflow encountered in double_scalars
  lambda x: np.log(p[0] * np.exp(-p[1] * x) + p[2]))
/home/mcf/anaconda3/lib/python3.5/site-packages/scipy/optimize/optimize.py:663: RuntimeWarning: invalid value encountered in double_scalars
  grad[k] = (f(*((xk + d,) + args)) - f0) / d[k]
Warning: Desired error not necessarily achieved due to precision loss.
         Current function value: -34.366008
         Iterations: 1
         Function evaluations: 552
         Gradient evaluations: 108


      fun: -34.36600756246744
 hess_inv: array([[ 0.99953426,  0.00164311, -0.04118523],
       [ 0.00164311,  0.99432624,  0.12712279],
       [-0.04118523,  0.12712279,  0.04248063]])
      jac: array([ -3.61366653,  10.76308966, -26.28632689])
  message: 'Desired error not necessarily achieved due to precision loss.'
     nfev: 552
      nit: 1
     njev: 108
   status: 2
  success: False
        x: array([3.07825766, 1.2641401 , 1.4789514 ])

### Using curve_fit

fit values: [2.55424137 1.35192223 0.4745096 ]
covariance matrix:
[[ 0.01588964  0.00681668 -0.00076153]
 [ 0.00681668  0.02019715  0.00541932]
 [-0.00076153  0.00541932  0.0028263 ]]
uncertainties: [0.12605411 0.14211667 0.05316297]

Мое наивное предположение, что в качестве параметра cв функции не допускает, чтобы NLL-функции функции были приятно представлены в некоторой форме, где нет экспоненты, свидетельствующей о том, что это приводит к сбою минимизации, поскольку исследуется некоторая область фазового пространства, которая вызывает inf s.

Если это предположение верно (?), То как вообще защититься от подобных проблем, учитывая, что эта игрушечная функция очень проста?Если мое предположение неверно, то что мне делать вместо этого?

...