Expectation of PERT RV is taking forever

[ritesh99rakesh]

I have the following implementation of PERT distribution (see PR (#16970) and #16970 (comment))

class PERTDistribution(SingleContinuousDistribution):
    _argnames = ('a', 'b', 'c')

    @property
    def set(self):
        return Interval(self.a, self.c)

    @staticmethod
    def check(a, b, c):
        _value_check((b > a, b < c),
            "Parameter b must be in range (%s, %s). b = %s"%(a, c, b))
        _value_check((a < c),
            "Parameter a must be less than Parameter c. a = %s, c = %s"%(a, c))

    def pdf(self, x):
        a, b, c = self.a, self.b, self.c
        alpha = (4*b + c - 5*a)/(c - a)
        beta = (5*c - a - 4*b)/(c - a)
        num = (x - a)**(alpha - 1)*(c - x)**(beta - 1)
        den = beta_fn(alpha, beta)*(c - a)**5
        return num/den

    # def expectation(self, expr, var, **kwargs):
    #     a, b, c = self.args
    #     return (a + 4*b + c)/6

But the following takes forever

>>> from sympy.stats import *
>>> froms sympy import symbols
>>> a, b, c = symbols('a b c')
>>> X = PERT('x', a, b, c)
>>> E(X) # takes forever without explicit definition of method `expectation`

Reason: Might be an issue with integrals (see #16970 (comment))

[jksuom]

It returns immediately on my system:

Python 2.7.6 (default, Nov 13 2018, 12:45:42) 
[GCC 4.8.4] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> from sympy.stats import *
>>> from sympy import symbols
>>> a, b, c = symbols('a, b, c')
>>> X = PERT('x', a, b, c)
>>> E(X)
a/6 + 2*b/3 + c/6

[ritesh99rakesh]

It returns immediately on my system:

Python 2.7.6 (default, Nov 13 2018, 12:45:42) 
[GCC 4.8.4] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> from sympy.stats import *
>>> from sympy import symbols
>>> a, b, c = symbols('a, b, c')
>>> X = PERT('x', a, b, c)
>>> E(X)
a/6 + 2*b/3 + c/6

@jksuom I have implemented explicit function for expectation in the above code. If you comment out the function as I did in the issue, it will take forever.

[smichr]

If the expectation is what you have commented out, why are you commenting it out?

[jksuom]

It seems that the incomplete beta integral is not implemented. (Long time is spent in meijerint and heurisch.)

[ritesh99rakesh]

@smichr In continuous distributions in sympy/stats/crv_types.py, the expectation is generally not implemented explicitly as it is just calculated by integrating over. Since this (implicit calculation) was taking a very long time, I implemented the function explicitly. But as suggested by @sidhantnagpal and I also believe that it might be an issue with integration method (see this #16970 (comment)). So I opened this issue.

Yes, @jksuom incomplete beta integral is not implemented and I am working on it. Should be able to send PR soon.

[oscarbenjamin]

Here is the integral:

from sympy import *
from sympy import beta as beta_fn
x, a, b, c = symbols('x, a, b, c')
alpha = (4*b + c - 5*a)/(c - a)
beta = (5*c - a - 4*b)/(c - a)
num = (x - a)**(alpha - 1)*(c - x)**(beta - 1)
den = beta_fn(alpha, beta)*(c - a)**5
integrand = num/den

pprint(integrand)
pprint(integrand.integrate((x, a, c))) # slow

The integrand looks like

             -5⋅a + 4⋅b + c             -a - 4⋅b + 5⋅c
        -1 + ──────────────        -1 + ──────────────
                 -a + c                     -a + c    
(-a + x)                   ⋅(c - x)                   
──────────────────────────────────────────────────────
             5  ⎛-5⋅a + 4⋅b + c  -a - 4⋅b + 5⋅c⎞      
     (-a + c) ⋅Β⎜──────────────, ──────────────⎟      
                ⎝    -a + c          -a + c    ⎠

I've put the integration as going from a to c. Does expectation know that the pdf is zero outside of [a, c] or will it be integrating from -oo to oo?