pythonscipymathematical-optimizationmlehessian-matrix

Retrieve optimization results from MLE by scipy.stats.fit()?

I am trying to estimate different distributions' parameters via scipy.stats."some_dist".fit(), and I'm having extreme difficulty retrieving any information on the optimization process being used. Specifically, I am looking for the Hessian, which most implemented algorithms output, as documented here.

For all its optimizations, SciPy returns an object called OptimizeResult, which contains useful info such as the Hessian, so how would one get this after .fit()'ing some data, which calls the same optimizations?

Solution

  • Looks like a bit of source-diving is needed; fortunately, it does not require copying a lot of source code. This is how basic fit works: 

    from scipy.stats import cauchy
data = [0, 3, 4, 4, 5, 9]
res = cauchy.fit(data)   # (3.9798237305661255, 0.9205374643383732)

    and this is how it is modified to return OptimizeResult:

    from scipy.optimize import minimize
args = cauchy._fitstart(data)
x0, func, restore, args = cauchy._reduce_func(args, {})
res = minimize(func, x0, args=(data,), method='BFGS')

    Now res is 

          fun: 14.337039523098689
 hess_inv: array([[ 0.23321703, -0.0117229 ],
       [-0.0117229 ,  0.36807373]])
      jac: array([ 0.0000000e+00, -1.1920929e-07])
  message: 'Optimization terminated successfully.'
     nfev: 32
      nit: 5
     njev: 8
   status: 0
  success: True
        x: array([3.9798262 , 0.92055376])

    where you may recognize the parameters as being res.x. The function being minimized is "penalized NNLF" (nonnegative likelihood function).

    By the way,

    For all its optimizations, SciPy returns an object called OptimizeResult

    is an over-generalization. This is true for minimize method. The scipy.stats.fit uses fmin by default, which returns no such thing.

    So, if want identical output to fmin, that can be arranged, but there'll be no extra information there.

    from scipy.optimize import fmin
args = cauchy._fitstart(data)
x0, func, restore, args = cauchy._reduce_func(args, {})
res = fmin(func, x0, args=(data,))
print(tuple(res))   #(3.9798237305661255, 0.9205374643383732)

    Using minimize with method='Nelder-Mead' has essentially the same effect. You do get a bit of extra information, but given that the method is simplex-based (no derivatives being computed), this information is of limited use. 

    res = minimize(func, x0, args=(data,), method='Nelder-Mead')
print(res)
print(tuple(res.x))

    prints

     final_simplex: (array([[3.97982373, 0.92053746],
       [3.97983057, 0.92060317],
       [3.97977536, 0.92059568]]), array([14.33703952, 14.33703953, 14.33703953]))
           fun: 14.337039523477827
       message: 'Optimization terminated successfully.'
          nfev: 58
           nit: 31
        status: 0
       success: True
             x: array([3.97982373, 0.92053746])
(3.9798237305661255, 0.9205374643383732)