使用 Python Sympy 求解方程组

如何解决使用 Python Sympy 求解方程组

我正在尝试使用 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)