时间戳

首页 > 解决方案 > 用 Python Sympy 求解方程组

python - 用 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, 0, 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), (i, 0, 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), (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)

这个最终命令返回“([], 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

谢谢！

标签： pythonsympysolverderivativeequation-solving

解决方案

基本上你想找到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

您正在尝试将数据替换为 forx[i]然后求解，但这并不是 sympy 的真正用途。fsolve如果您只是在寻找数值解决方案，最好使用 numpy 中的 etc。

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，所以我猜不是你想要的。

现在您可以将您的值替换为 forx[i]或将上面的公式转换为代码或使用lambdify或您想做的任何事情，例如：

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

推荐阅读