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 でこの積分を解く方法について考えている人はいますか?」のような質問にまとめることができます。

また、報告されたエラーの性質を考えると、あなたが遭遇した端末の問題は浮動小数点数の多項式にあったように思えます。5.5SymPy に次のように入力すると、あまり数学的ではない Python の浮動小数点数が生成されます。ただし、入力がきれいな整数であっても、SymPy はこの問題を解決できないため、別の問題があります。

また、一般的な注意点として、分析的に簡単に解決できない統計の問題を見つけるのは簡単です。SymPy が分析的に解決できない統計上の問題を見つけるのはさらに簡単です。ここでの私の経験則は、通常、学部の教科書に出てくる問題はすべて解決できるということです。

于 2014-04-01T14:52:34.453 に答える