Mixture of inverse gamma distributions in sympy.stats

Question

Following up on a previous post, I'd like to do the following to take a weighted mixture of inverse gamma distributions using 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]

Unfortunately this errors at D3 = density(G3), with an error that concludes in

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.

Is there a way round this?

Ta.

Answer 1

Sympy.stats is producing an integral that looks like the following:

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                                                                            

So one could reduce your question into a question like, "Does anyone have any thoughts on how SymPy could solve this integral?"

Also, given the nature of the reported error it sounds like maybe the terminal issue that you ran into was with polynomials over floating point numbers. Typing things like 5.5 into SymPy produces Python floats which aren't very mathy. SymPy can't solve this problem even when the inputs are clean integers though, so there are other issues.

Also just a general note, it's easy to find problems in statistics that aren't easily solvable analytically. It's even easier to find problems in statistics that SymPy can't solve analytically. My rule of thumb here is that we can usually solve any problem that you would find in an undergraduate textbook.