如何解决不能在 sympy 中整合简单的正态分布，这取决于均值和偏差常数
所以...我可以sympy.integrate具有均值和标准差的正态分布：
(10.1,0.333333333),# Works fine
但不是：
(8.655555555555557,0.5212875796916135),# Fails
种类感觉应该不会有太大的不同。 这是怎么回事？
完整示例：
import numpy as np
%matplotlib inline
import sympy
from sympy import symbols
from sympy import plot
def normal(x,mean,sigma):
z = (x - mean) / sigma
return (1 / (sigma * sympy.sqrt(2 * sympy.pi))) * sympy.exp(-(z * z) / 2)
for μ,σ in [
(10.1,# Works fine
(8.655555555555557,# Fails
]:
x = symbols("x")
print(f"μ={μ},σ={σ}")
distrib = normal(x=x,mean=μ,sigma=σ)
# distrib = sympy.simplify(distrib) # Doesn't help
distrib_cum = sympy.integrate(distrib,x)
# distrib_cum = sympy.simplify(distrib_cum) # Doesn't help
print(distrib_cum)
plot(distrib_cum)
NameError Traceback (most recent call last)
<ipython-input-18-e85bdeabeaa6> in <module>
22 # distrib_cum = sympy.simplify(distrib_cum) # Doesn't help
23 print(distrib_cum)
---> 24 plot(distrib_cum)
25
26
~/venv/p37_default/lib/python3.7/site-packages/sympy/plotting/plot.py in plot(show,*args,**kwargs)
1738 plots = Plot(*series,**kwargs)
1739 if show:
-> 1740 plots.show()
1741 return plots
1742
~/venv/p37_default/lib/python3.7/site-packages/sympy/plotting/plot.py in show(self)
220 self._backend.close()
221 self._backend = self.backend(self)
--> 222 self._backend.show()
223
224 def save(self,path):
~/venv/p37_default/lib/python3.7/site-packages/sympy/plotting/plot.py in show(self)
1414
1415 def show(self):
-> 1416 self.process_series()
1417 #Todo after fixing https://github.com/ipython/ipython/issues/1255
1418 # you can uncomment the next line and remove the pyplot.show() call
~/venv/p37_default/lib/python3.7/site-packages/sympy/plotting/plot.py in process_series(self)
1411 if isinstance(self.parent,PlotGrid):
1412 parent = self.parent.args[i]
-> 1413 self._process_series(series,ax,parent)
1414
1415 def show(self):
~/venv/p37_default/lib/python3.7/site-packages/sympy/plotting/plot.py in _process_series(self,series,parent)
1239 # Create the collections
1240 if s.is_2Dline:
-> 1241 collection = self.LineCollection(s.get_segments())
1242 ax.add_collection(collection)
1243 elif s.is_contour:
~/venv/p37_default/lib/python3.7/site-packages/sympy/plotting/plot.py in get_segments(self)
704 list_segments.append([p,q])
705
--> 706 f_start = f(self.start)
707 f_end = f(self.end)
708 sample(np.array([self.start,f_start]),~/venv/p37_default/lib/python3.7/site-packages/sympy/plotting/experimental_lambdify.py in __call__(self,args)
173 try:
174 #The result can be sympy.Float. Hence wrap it with complex type.
--> 175 result = complex(self.lambda_func(args))
176 if abs(result.imag) > 1e-7 * abs(result):
177 return None
~/venv/p37_default/lib/python3.7/site-packages/sympy/plotting/experimental_lambdify.py in __call__(self,**kwargs)
270
271 def __call__(self,**kwargs):
--> 272 return self.lambda_func(*args,**kwargs)
273
274
<string> in <lambda>(x0)
NameError: name 'Integral' is not defined
第一个问题：
μ=10.1,σ=0.333333333
0.353553390593274*sqrt(2)*erf(2.12132034568096*x - 21.4253354913777)
工作正常。
第二个：
μ=8.655555555555557,σ=0.5212875796916135
1.30203374788369e-60*sqrt(2)*Integral(exp(31.8522556903367*x)*exp(-1.83998909636091*x**2),x)/sqrt(pi)
绘制时出现“NameError: name 'Integral' is not defined' 失败。
（而且我有点怀疑 1.30203374788369e-60* .... 的结局是否会很好）。
截至今天的最新 pip 安装失败（python 3.7，sympy 1.18）。
我的错误还是同情的局限性？
解决方法
这是一个有效的关闭案例：
In [125]: normal(x,8.6,0.5)
Out[125]:
2
-147.92⋅(0.116279069767442⋅x - 1)
1.0⋅√2⋅ℯ
──────────────────────────────────────────
√π
In [126]: integrate(_,x)
Out[126]: 0.353553390593274⋅√2⋅erf(1.41421356237309⋅x - 12.1622366364086)
还有一个没有：
In [127]: normal(x,8.655555555555557,0.5212875796916135)
Out[127]:
2
-137.849484348735⋅(0.115532734274711⋅x - 1)
0.95916346270094⋅√2⋅ℯ
─────────────────────────────────────────────────────────────────
√π
In [128]: integrate(_,x)
Out[128]:
⌠
⎮ 2
⎮ 31.8522556903367⋅x -1.83998909636091⋅x
1.30203374788369e-60⋅√2⋅⎮ ℯ ⋅ℯ dx
⌡
──────────────────────────────────────────────────────────────────────
√π
我不知道为什么会有不同。有时在 sympy 情况下，我们需要向符号添加约束，例如 real 或 positive。看起来它无法将积分减少到 erf。我想知道是否存在 erf 定义有效的值范围。我已经很多年没接触过它了。
正如您所注意到的，前导 1.30203374788369e-60 是可疑的，尤其是与工作案例中的 1 相比。
从 sigma 中去掉几个数字就可以了：0.52128757969161
TLDR：在可以避免的情况下不要使用浮点数。 Python 的有限精度与 SymPy 的无限精度算法不能很好地配合。
这与0.1 + 0.2 = 0.30000000000000004的故事非常相似。
对于这两种情况，SymPy 最终都会尝试使用 Meijer G 算法来评估积分。在它自动删除一个系数后，它会尝试对表达式进行积分。
因此，对于第一个示例，它尝试将 exp(-459.04500091809*(0.099009900990099*x - 1)**2) 与替换 y = x + 10.1000000000000 集成。这导致了 exp(-4.500000009*x**2)，Meijer 可以解决这个问题。
但是，对于第二个示例，它尝试将 exp(-137.849484348735*(0.115532734274711*x - 1)**2) 与替换 y = x + 8.65555555555556 集成。这导致了 exp(-1.83998909636091*(x - 9.60959706869997e-16)**2)，Meijer 无法求解，因为它没有 exp(-x^2) 的典型形式。
有3种方法可以解决问题：
- 使用具有无限精度的东西，例如有理数。根据计算类型和试验值的数量，此或选项 2 可能是最有效的：
import sympy
from sympy import symbols
from sympy import plot
def normal(x,mean,sigma):
z = (x - mean) / sigma
return (1 / (sigma * sympy.sqrt(2 * sympy.pi))) * sympy.exp(-(z * z) / 2)
x = symbols("x")
for μ,σ in [
(sympy.Rational("10.1"),sympy.Rational(1,3)),# Works fine
(sympy.Rational("8.655555555555557").limit_denominator(10**10),sympy.Rational("0.5212875796916135").limit_denominator(10**10)),# Fails
]:
print(f"μ={μ},σ={σ}")
distrib = normal(x,μ,σ)
distrib_cum = sympy.integrate(distrib,x)
print(distrib_cum)
plot(distrib_cum)
- 计算一次，然后替换测试值。在某些情况下，您可能需要将某些符号指定为实数或正数。
import sympy
from sympy import symbols
from sympy import plot
def normal(x,sigma):
z = (x - mean) / sigma
return (1 / (sigma * sympy.sqrt(2 * sympy.pi))) * sympy.exp(-(z * z) / 2)
x,mu,sigma = symbols("x,sigma")
distrib = normal(x,sigma)
distrib_cum = sympy.integrate(distrib,x)
for μ,σ in [
(10.1,0.333333333),# Works fine
(8.655555555555557,0.5212875796916135),σ={σ}")
print(distrib_cum.subs({mu: μ,sigma: σ}))
plot(distrib_cum.subs({mu: μ,sigma: σ}))
- 也许 SymPy 已经为您计算过了。当然，这是cdf而不是pdf的反导数。常量之间存在细微差别。
from sympy import symbols
from sympy import plot
from sympy.stats import Normal,cdf
x = symbols('x')
for μ,σ={σ}")
distrib_cum = cdf(Normal('X',σ))(x)
print(distrib_cum)
plot(distrib_cum)
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。