scipy - Evaluate Derivative of Negative Log Likelihood Python

Question

I am trying to evaluate the derivative of the negative log likelihood functionin python. I am using sympy to compute the derivative however, I receive an error when I try to evaluate it. The example code below attempts to compute this for the lognormal function.

import numpy as np
import sympy as sym
from sympy import Product, Function, oo, IndexedBase, diff, Eq, symbols, log, exp
from scipy.stats import lognorm

np.random.seed(seed=111)
test = lognorm.rvs(s=1,loc=2,scale=1,size=1000)

x = IndexedBase('x')
i = symbols('i', positive=True)
n = symbols('n', positive=True)
mx = symbols('mx', positive=True)
sx = symbols('sx', positive=True)

pdf = 1 / (x[i] * sx * sqrt(2 * np.pi)) * exp(-0.5 * ((log(x[i]-mx))**2/(sx)**2))
Log_LL = -log(Product(pdf, (i, 1, n)))

deriv = diff(Log_LL, mx) 
fx = lambdify([x,mx,sx,n],deriv)
fx(test,2,1,len(test))

When I evaluate this formula, I receive the following error:

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-139-452f00dee1f1> in <module>
      1 fx = lambdify([x,mx,sx,n],deriv)
----> 2 fx(test,2,1,len(test))

<lambdifygenerated-14> in _lambdifygenerated(Dummy_39, mx, sx, n)
      3   # Derivative
      4   # Product
----> 5 -Derivative(Product(0.398942280401433*exp(-0.5*log(-mx + Dummy_39[i])**2/sx**2)/(sx*Dummy_39[i]), (i, 1, n)), mx)/Product(0.398942280401433*exp(-0.5*log(-mx + Dummy_39[i])**2/sx**2)/(sx*Dummy_39[i]), (i, 1, n)))

NameError: name 'Derivative' is not defined

I realize that the expression "deriv" contains the derivative operator, but I am simply computing the derivative of a single variable, so I believe sympy should be able to handle this.

I am running sympy 1.7.1, numpy 1.19.2 and scipy 1.5.2

Thank you!

question from:https://stackoverflow.com/questions/66066327/evaluate-derivative-of-negative-log-likelihood-python

Answer 1

First off, the symbolic sympy and the numeric numpy/scipy don't mix well. It is usually recommended separating them. Also, sympy doesn't like floats as these prevent exact symbolic calculations. Therefore, instead of np.pi, sympy's symbolic pi is preferred. Also, 0.5 should be replaced by a sympy fraction, for example S.Half. (See also sympy's gotchas.)

Sympy's derivative doesn't seem to be able to cope with the Product. We can try to replace the log of the product by a sum of the logs. expand_log(..., force=True) can help with that conversion (force=True when sympy isn't sure that the expression is certain to be positive, presumably the x[i] could be complex).

When converting to numpy, numpy doesn't like the indexing starting from 1. This can be solved with indexing starting from 0.

from sympy import Product, Sum, IndexedBase, diff, symbols, log, exp, lambdify, sqrt, pi, S, expand_log

x = IndexedBase('x')
i = symbols('i', positive=True)
n = symbols('n', positive=True)
mx = symbols('mx', positive=True)
sx = symbols('sx', positive=True)

pdf = 1 / (x[i] * sx * sqrt(2 * pi)) * exp(-S.Half * ((log(x[i] - mx)) ** 2 / (sx) ** 2))
# Log_LL = -Sum(log(pdf), (i, 0, n - 1))
Log_LL = -log(Product(pdf, (i, 0, n-1)))
Log_LL = expand_log(Log_LL, force=True)

deriv = diff(Log_LL, mx)
fx = lambdify([x, mx, sx, n], deriv)

from scipy.stats import lognorm
import numpy as np

np.random.seed(seed=111)
test = lognorm.rvs(s=1, loc=2, scale=1, size=1000)
fx(test, 2, 1, len(test))

This way, sympy manages to calculate the derivative symbolically:

 n - 1                 
  ____                 
  ╲                    
   ╲   log(-mx + x[i]) 
    ╲  ────────────────
-   ╱    2             
   ╱   sx ?(-mx + x[i])
  ╱                    
  ￣￣￣￣                 
 i = 0