用 Python Sympy 求解方程组

Question

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

这条最终命令returns“([], 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

谢谢！

Answer 1

基本上你想找到 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], [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                         ⎠

另一个解决方案 s1 是否定的 sigma 所以我想这不是你想要的。

现在您可以将您的值替换为 x[i] 或将上面的公式转换为代码或使用 lambdify 或您想要执行的任何操作，例如：

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

用 Python Sympy 求解方程组

Solve System of Equations with Python Sympy

python

solver

sympy

derivative

equation-solving