问题描述
我正在尝试使用 Sympy 在 Python 中求解由两个方程组成的系统。这比标准问题有点棘手，因为它包含两个方程的求和，其中两个方程都是通过对对数正态 pdf 的负对数似然求导来找到的。这是我的代码：
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,n-1)))
Log_LL22 = expand_log(Log_LL2,force=True)
pprint(Log_LL22)
返回：
-Sum(-log(sigma) - log(x[i]) - log(pi)/2 - log(2)/2 - (-mu + log(x[i]))**2/(2*sigma**2),n - 1))
df_dmu = diff(Log_LL22,mu)
df_dsigma = diff(Log_LL22,sigma)
pprint(df_dmu )
pprint(df_dsigma )
返回：
-Sum(-(2*mu - 2*log(x[i]))/(2*sigma**2),n - 1))
-Sum(-1/sigma + (-mu + log(x[i]))**2/sigma**3,n - 1))
solve([df_dmu.subs([(n,len(test['y'])),(x,test['y'])]),df_dsigma.subs([(n,test['y'])])],mu,sigma,set=True)
最后一条命令返回“([],set())”。我不确定如何实现这个方程组，同时告诉求解器替换 x_i 和 n 以求解 mu 和 sigma。如果可能的话，我也很乐意不插入 x_i 和 n 并根据 x_i 和 n 接收答案。我知道这些参数可以用 scipy 的拟合函数解决，但是，当我计算负对数似然的 Hessian 并插入 scipy 拟合参数时，结果是 scipy 拟合参数和手动计算参数之间的数量级差异。
我正在运行 sympy 1.7.1、numpy 1.19.2 和 scipy 1.5.2
谢谢！
解决方法
基本上你想找到 mu 和 sigma 使这些表达式为零：
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
您试图用 x[i] 中的数据代替，然后求解，但这并不是 sympy 的真正用途。如果您只是在寻找数值解，最好使用 numpy 中的 fsolve 等。
sympy 是为了找到解决方案的通用公式。所以我们只想解决上面的 mu 和 sigma。不幸的是，solve 不明白如何处理这样的求和，所以它放弃了：
In [36]: solve([df_dmu,df_dsigma],[mu,sigma])
Out[36]: []
我们可以通过操纵方程从求和中提取感兴趣的符号来解决问题：
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
───────────────────────────────────────────────────────────────
σ
现在solve给出一般解决方案：
In [54]: s1,s2 = solve([eq_mu,eq_sigma],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 ⎠
另一个解决方案 s1 是针对负 sigma 所以我猜不是你想要的。
现在您可以将您的值替换为 x[i] 或将上面的公式转换为代码或使用 lambdify 或您想做的任何事情，例如：
In [57]: lambdify((n,x),s2)(3,[5,6,7])
Out[57]: (1.7823691769058227,0.1375246031240532)