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!
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)