error when taking mixture of inverse gamma distributions

[John Griffiths]

Does anyone have any thoughts on how to solve this problem: 

When I try to take a weighted mixture of inverse gamma distributions using sympy.stats I get the following error

%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]

...which errors at  D3 = density(G3). The error includes

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.

(also get this when I take mixture of inverse Gamma with Normal and Uniform distributions)

Should this be possible?

Cheers. 

(p.s. apologies for redundancy with recent SO post)

[Matthew Rocklin]

First, thanks for posting the question on stack overflow, I think that we should use that more.  Unfortunately, few of us (maybe only Aaron?) actually checks SO, again, I think that we should use it more.

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?"

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/4cb36d1b-6827-4123-8efe-b1462d822e6b%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Aaron Meurer]

Yes, you can go to http://stackexchange.com/filters/ and set up a filter for the SymPy tag, and it will email you new questions. Lately there are new questions almost every single day, so it would definitely help if others were answering there as well.

The questions on stackoverflow tend to be high quality and often have rather interesting use-cases (there are doozies, of course, but the worst questions are closed automatically).

Aaron Meurer

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAJ8oX-HEDx4fraBRj4HARPaCEJKXKaTP3zjS_dbeqd3WoeAAWA%40mail.gmail.com.

[Matthew Rocklin]

Done.  You should announce that tip to the list.  There were indeed some good questions on there.  SO is a good way to keep our collective ear to the ground.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6JUTWn5vz7oih7%2BY18j6jNbvxZD4uEy0L9zTtp4-f3www%40mail.gmail.com.

[Aaron Meurer]

This is the list...

Aaron Meurer

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAJ8oX-EUvvPMJhPMb6FqoyuxY9MS2OLGavnHavgH0dQLBaBjpQ%40mail.gmail.com.

[Matthew Rocklin]

It's to anyone who reads a message with title about inverse gamma distributions.  Which I suspect is just John, because he asked the question, myself, because I curate stats, and you, because you read everything :)

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6Kb83%2B9nBkPJpCHFEBMxzNHgtPjo%3DKFpsTW6yh3P-hc6w%40mail.gmail.com.