Solve System of Equations with Python Sympy

Question

I am trying to solve a system of two equations in Python using Sympy. This is bit trickier than a standard problem because it includes a summation in both equations, where both equations are found by taking the derivative of the negative loglikelihood of a lognormal pdf. Here is my code:

import numpy as np
from pprint import pprint
import sympy as sym
from sympy.solvers import solve
from sympy import Product, Function, oo, IndexedBase, diff, Eq, symbols, log, exp, pi, S, expand_log
from scipy.stats import lognorm
import scipy

np.random.seed(seed=111)
test = pd.DataFrame(data=lognorm.rvs(s=1,loc=2,scale=1,size=1000),columns=['y'])

x = IndexedBase('x')
i = symbols('i', positive=True)
n = symbols('n', positive=True)
mu = symbols('mu', positive=True)
sigma = symbols('sigma', positive=True)

pdf2 = 1 / (x[i] * sigma * sqrt(2 * pi)) * exp(-S.Half * ((log(x[i])-mu)/(sigma))**2)
Log_LL2 = -log(Product(pdf2, (i, 0, n-1)))
Log_LL22 = expand_log(Log_LL2, force=True)
pprint(Log_LL22)

Returns:

-Sum(-log(sigma) - log(x[i]) - log(pi)/2 - log(2)/2 - (-mu + log(x[i]))**2/(2*sigma**2), (i, 0, n - 1))

df_dmu = diff(Log_LL22, mu)
df_dsigma = diff(Log_LL22, sigma)
pprint(df_dmu )
pprint(df_dsigma )

Returns:

-Sum(-(2*mu - 2*log(x[i]))/(2*sigma**2), (i, 0, n - 1))
-Sum(-1/sigma + (-mu + log(x[i]))**2/sigma**3, (i, 0, n - 1))

solve([df_dmu.subs([(n,len(test['y'])),(x,test['y'])]),df_dsigma.subs([(n,len(test['y'])),(x,test['y'])])],mu,sigma,set=True)

This final command returns "([], set())". I am not sure how to implement this system of equations while telling the solver to substitute for x_i and n in order to solve for mu and sigma. I would also be happy not to plug in x_i and n and receive the answer in terms of x_i and n if this was possible. I understand these parameters can be solved for with scipy's fit function, however, when I compute the Hessian of the negative loglikelihood and plug in the scipy fitted parameters, the result is orders of magnitude difference between the scipy fitted parameters and the manually computed parameters.

I am running sympy 1.7.1, numpy 1.19.2 and scipy 1.5.2

Thank you!

Answer 1

Basically you want to find mu and sigma to make these expressions zero:

In [47]: df_dmu
Out[47]: 
 n - 1                      
  ____                      
  ╲                         
   ╲   -(2⋅μ - 2⋅log(x[i])) 
    ╲  ─────────────────────
-   ╱              2        
   ╱            2⋅σ         
  ╱                         
  ‾‾‾‾                      
 i = 0                      

In [48]: df_dsigma
Out[48]: 
 n - 1                          
 _____                          
 ╲                              
  ╲                             
   ╲   ⎛                      2⎞
    ╲  ⎜  1   (-μ + log(x[i])) ⎟
-   ╱  ⎜- ─ + ─────────────────⎟
   ╱   ⎜  σ            3       ⎟
  ╱    ⎝              σ        ⎠
 ╱                              
 ‾‾‾‾‾                          
 i = 0    

You're trying to substitute the data in for x[i] and then solve but that's not really what sympy is for. You'd be better of using fsolve etc from numpy if you're just looking for a numerical solution.

What sympy is for is for finding the general formula for the solution. So we just want to solve the above for mu and sigma. Unfortunately solve doesn't understand how to work with summations like this so it gives up:

In [36]: solve([df_dmu, df_dsigma], [mu, sigma])
Out[36]: []

We can give solve a hand though by manipulating the equations to extract the symbols of interest from the summations:

In [49]: eq_mu = factor_terms(expand(df_dmu))

In [50]: eq_sigma = factor_terms(expand(df_dsigma))

In [51]: eq_mu
Out[51]: 
  n - 1     n - 1          
   ___       ___           
   ╲         ╲             
    ╲         ╲            
μ⋅  ╱   1 -   ╱   log(x[i])
   ╱         ╱             
   ‾‾‾       ‾‾‾           
  i = 0     i = 0          
───────────────────────────
              2            
             σ             

In [52]: eq_sigma
Out[52]: 
     n - 1         n - 1                       n - 1           
      ___           ___                         ___            
      ╲             ╲                           ╲              
   2   ╲             ╲                           ╲      2      
  μ ⋅  ╱   1   2⋅μ⋅  ╱   log(x[i])   n - 1       ╱   log (x[i])
      ╱             ╱                 ___       ╱              
      ‾‾‾           ‾‾‾               ╲         ‾‾‾            
     i = 0         i = 0               ╲       i = 0           
- ────────── + ─────────────────── +   ╱   1 - ────────────────
       2                 2            ╱                2       
      σ                 σ             ‾‾‾             σ        
                                     i = 0                     
───────────────────────────────────────────────────────────────
                               σ   

Now solve gives the general solution:

In [54]: s1, s2 = solve([eq_mu, eq_sigma], [mu, sigma])

In [55]: s2
Out[55]: 
⎛                           _________________________________________⎞
⎜                          ╱                                       2 ⎟
⎜n - 1                    ╱    n - 1              ⎛n - 1          ⎞  ⎟
⎜ ___                    ╱      ___               ⎜ ___           ⎟  ⎟
⎜ ╲                     ╱       ╲                 ⎜ ╲             ⎟  ⎟
⎜  ╲                   ╱         ╲      2         ⎜  ╲            ⎟  ⎟
⎜  ╱   log(x[i])      ╱      n⋅  ╱   log (x[i]) - ⎜  ╱   log(x[i])⎟  ⎟
⎜ ╱                  ╱          ╱                 ⎜ ╱             ⎟  ⎟
⎜ ‾‾‾               ╱           ‾‾‾               ⎜ ‾‾‾           ⎟  ⎟
⎜i = 0            ╲╱           i = 0              ⎝i = 0          ⎠  ⎟
⎜───────────────, ───────────────────────────────────────────────────⎟
⎝       n                                  n                         ⎠

The other solution s1 is for negative sigma so I guess not what you want.

Now you can substitute your values in for x[i] or translate the formula above into code or use lambdify or whatever it is you want to do e.g.:

In [57]: lambdify((n, x), s2)(3, [5, 6, 7])
Out[57]: (1.7823691769058227, 0.1375246031240532)