0

继上一篇文章之后，我想执行以下操作以使用 sympy.stats 对逆伽马分布进行加权混合：

%matplotlib inline
from matplotlib import pyplot as plt
from sympy.stats import GammaInverse, density
import numpy as np

f1 = 0.7; f2 = 1-f1
G1 = GammaInverse("G1", 5, 120/(5.5*2.5E-7))
G2 = GammaInverse("G2", 4, 120/(5.5*1.5E-7))
G3 = f1*G1 + f2*G2
D1 = density(G1);  
D2 = density(G2);  
D3 = density(G3);
v1 = [D1(i).evalf() for i in u]
v2 = [D2(i).evalf() for i in u]
v3 = [D3(i).evalf() for i in u]

不幸的是，这个错误在D3 = density(G3)，错误结束于

PolynomialDivisionFailed: couldn't reduce degree in a polynomial 
division algorithm when dividing [231761.370742578/(0.0011381138741823*G2**2 -
 0.007587425827882*G2*_z + 0.0126457097131367*_z**2), 0.0]
by [263.770831541635/263.770831541635, 0.0]. 
This can happen when it's not possible to detect zero in the coefficient domain. 
The domain of computation is RR(G2,_t0,_z). Zero detection is guaranteed in this
coefficient domain. This may indicate a bug in SymPy or the domain is user defined
and doesn't implement zero detection properly.

有办法解决这个问题吗？

塔。

4

1 回答 1

0

Sympy.stats 正在生成如下所示的积分：

In [1]: from sympy.stats import *
In [2]: a, b, c, d, = symbols('a b c d', real=True, positive=True)
In [3]: G1 = GammaInverse("G1", a, b)
In [4]: G2 = GammaInverse("G2", c, d)
In [5]: G3 = S(7)/10*G1 + S(3)/10*G2
In [7]: density(G3, evaluate=False)(x)
Out[7]: 
∞                                                                            
⌠                                                                            
⎮                  ∞                                                         
⎮                  ⌠                                                         
⎮                  ⎮              -b                                         
⎮                  ⎮              ───                                        
⎮              -d  ⎮   -a - 1  a   G₁           ⎛7⋅G₁   3⋅G₂    ⎞            
⎮              ─── ⎮ G₁      ⋅b ⋅ℯ   ⋅DiracDelta⎜──── + ──── - x⎟            
⎮   -c - 1  c   G₂ ⎮                            ⎝ 10     10     ⎠            
⎮ G₂      ⋅d ⋅ℯ   ⋅⎮ ──────────────────────────────────────────── d(G₁)      
⎮                  ⎮                     Γ(a)                                
⎮                  ⌡                                                         
⎮                  0                                                         
⎮ ───────────────────────────────────────────────────────────────────── d(G₂)
⎮                                  Γ(c)                                      
⌡                                                                            
0

因此，可以将您的问题简化为一个问题，例如“有人对 SymPy 如何解决这个积分有任何想法吗？”

此外，鉴于报告错误的性质，听起来您遇到的终端问题可能是浮点数上的多项式。在 SymPy 中输入类似的东西5.5会产生不是很数学的 Python 浮点数。即使输入是干净的整数，SymPy 也无法解决这个问题，所以还有其他问题。

也只是一般性说明，很容易发现统计中不易通过分析解决的问题。更容易找到 SymPy 无法解析解决的统计问题。我的经验法则是，我们通常可以解决您在本科教科书中发现的任何问题。

于 2014-04-01T14:52:34.453 回答