用Python SymPy求解含求和项的对数正态似然方程组空结果问题
SymPy求解带求和项的对数正态似然方程问题
场景说明
在Python环境下借助SymPy求解一组二元方程,两个方程均由对数正态概率密度函数的负对数似然求导得到,内部包含求和项,求解复杂度高于常规场景。
测试代码
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))
求导后得到的目标方程
对mu和sigma分别求导的代码如下:
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))
遇到的问题
- 执行如下求解命令时,最终返回空结果
([], set())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) - 不清楚如何在求解时告知求解器替换求和项中的x_i与n变量,以得到mu和sigma的数值解
- 希望能够不预先代入x_i和n的具体值,直接得到以x_i和n表示的解析形式解
- 已知可以通过scipy的fit函数求解上述参数,但计算负对数似然的Hessian矩阵并代入scipy拟合得到的参数时,结果显示scipy拟合参数与手动计算的参数存在数个数量级的差异
运行环境版本
- sympy 1.7.1
- numpy 1.19.2
- scipy 1.5.2
内容的提问来源于stack exchange,提问作者njalex22
相关产品推荐
相关产品推荐

