Repository: sympy/sympy

Issue:
Error with handmade Gamma Distribution (Easy Fix)
## Problem
The following was working in sympy version 1.5. It does not work in 1.6.

The following code returns `ValueError: The pdf on the given set is incorrect.`:
```python
from sympy.stats import *
from sympy import *

x, k, theta = symbols(r"x k theta", positive=True)
pdf = 1/(gamma(k)*theta**k)*x**(k-1)*exp(-x/theta)

X = ContinuousRV(x, pdf, set=Interval(0, oo))
```
However, this is exactly the pdf in the documentation.

Upon inspection, I found in [line 154](https://github.com/sympy/sympy/blob/8f6bc5947262754fe623d30682e47b7e01580bf7/sympy/stats/crv_types.py#L154) in `crv_types.py`, we have that `val = integrate(pdf(x), (x, set))` returns `theta*theta**(-k)*theta**(k - 1)`. I assume `integrate` does not simplify the expression since this value is 1.

## Proposed possible fixes:
- Replace `val = integrate(pdf(x), (x, set))` with `val = simplify(integrate(pdf(x), (x, set)))`. (This may be inefficient for complex integrals that turn out to be incorrect).
- Replace `_value_check(val == S.One, "The pdf on the given set is incorrect.")` with `_value_check((val - S.One) == S.Zero, "The pdf on the given set is incorrect.")`. (This might also use the `simplify` method behind the scenes. But if it does not, it will be more efficient).

`drv_types.py` may have the same issue.

Unless I have done something horribly wrong, should I submit a PR for this? (I'll also add unit tests for this).

Discussion:

Patch:
diff --git a/sympy/stats/crv.py b/sympy/stats/crv.py
index 16dda9c70f..fa568d2a29 100644
--- a/sympy/stats/crv.py
+++ b/sympy/stats/crv.py
@@ -333,7 +333,7 @@ def compute_characteristic_function(self, **kwargs):
         """
         x, t = symbols('x, t', real=True, cls=Dummy)
         pdf = self.pdf(x)
-        cf = integrate(exp(I*t*x)*pdf, (x, -oo, oo))
+        cf = integrate(exp(I*t*x)*pdf, (x, self.set))
         return Lambda(t, cf)
 
     def _characteristic_function(self, t):
@@ -355,7 +355,7 @@ def compute_moment_generating_function(self, **kwargs):
         """
         x, t = symbols('x, t', real=True, cls=Dummy)
         pdf = self.pdf(x)
-        mgf = integrate(exp(t * x) * pdf, (x, -oo, oo))
+        mgf = integrate(exp(t * x) * pdf, (x, self.set))
         return Lambda(t, mgf)
 
     def _moment_generating_function(self, t):
diff --git a/sympy/stats/crv_types.py b/sympy/stats/crv_types.py
index 958f43dd5f..12a1774aff 100644
--- a/sympy/stats/crv_types.py
+++ b/sympy/stats/crv_types.py
@@ -127,10 +127,11 @@
 def _(x):
     return any([is_random(i) for i in x])
 
-def rv(symbol, cls, args):
+def rv(symbol, cls, args, **kwargs):
     args = list(map(sympify, args))
     dist = cls(*args)
-    dist.check(*args)
+    if kwargs.pop('check', True):
+        dist.check(*args)
     pspace = SingleContinuousPSpace(symbol, dist)
     if any(is_random(arg) for arg in args):
         from sympy.stats.compound_rv import CompoundPSpace, CompoundDistribution
@@ -152,10 +153,10 @@ def set(self):
     def check(pdf, set):
         x = Dummy('x')
         val = integrate(pdf(x), (x, set))
-        _value_check(val == S.One, "The pdf on the given set is incorrect.")
+        _value_check(Eq(val, 1) != S.false, "The pdf on the given set is incorrect.")
 
 
-def ContinuousRV(symbol, density, set=Interval(-oo, oo)):
+def ContinuousRV(symbol, density, set=Interval(-oo, oo), **kwargs):
     """
     Create a Continuous Random Variable given the following:
 
@@ -168,6 +169,11 @@ def ContinuousRV(symbol, density, set=Interval(-oo, oo)):
         Represents probability density function.
     set : set/Interval
         Represents the region where the pdf is valid, by default is real line.
+    check : bool
+        If True, it will check whether the given density
+        integrates to 1 over the given set. If False, it
+        will not perform this check. Default is False.
+
 
     Returns
     =======
@@ -196,7 +202,9 @@ def ContinuousRV(symbol, density, set=Interval(-oo, oo)):
     """
     pdf = Piecewise((density, set.as_relational(symbol)), (0, True))
     pdf = Lambda(symbol, pdf)
-    return rv(symbol.name, ContinuousDistributionHandmade, (pdf, set))
+    # have a default of False while `rv` should have a default of True
+    kwargs['check'] = kwargs.pop('check', False)
+    return rv(symbol.name, ContinuousDistributionHandmade, (pdf, set), **kwargs)
 
 ########################################
 # Continuous Probability Distributions #
diff --git a/sympy/stats/drv_types.py b/sympy/stats/drv_types.py
index c0fb01946f..7018014e88 100644
--- a/sympy/stats/drv_types.py
+++ b/sympy/stats/drv_types.py
@@ -17,7 +17,7 @@
 
 from sympy import (Basic, factorial, exp, S, sympify, I, zeta, polylog, log, beta,
                    hyper, binomial, Piecewise, floor, besseli, sqrt, Sum, Dummy,
-                   Lambda)
+                   Lambda, Eq)
 from sympy.stats.drv import SingleDiscreteDistribution, SingleDiscretePSpace
 from sympy.stats.rv import _value_check, is_random
 
@@ -33,10 +33,11 @@
 ]
 
 
-def rv(symbol, cls, *args):
+def rv(symbol, cls, *args, **kwargs):
     args = list(map(sympify, args))
     dist = cls(*args)
-    dist.check(*args)
+    if kwargs.pop('check', True):
+        dist.check(*args)
     pspace = SingleDiscretePSpace(symbol, dist)
     if any(is_random(arg) for arg in args):
         from sympy.stats.compound_rv import CompoundPSpace, CompoundDistribution
@@ -58,9 +59,9 @@ def set(self):
     def check(pdf, set):
         x = Dummy('x')
         val = Sum(pdf(x), (x, set._inf, set._sup)).doit()
-        _value_check(val == S.One, "The pdf is incorrect on the given set.")
+        _value_check(Eq(val, 1) != S.false, "The pdf is incorrect on the given set.")
 
-def DiscreteRV(symbol, density, set=S.Integers):
+def DiscreteRV(symbol, density, set=S.Integers, **kwargs):
     """
     Create a Discrete Random Variable given the following:
 
@@ -73,6 +74,10 @@ def DiscreteRV(symbol, density, set=S.Integers):
         Represents probability density function.
     set : set
         Represents the region where the pdf is valid, by default is real line.
+    check : bool
+        If True, it will check whether the given density
+        integrates to 1 over the given set. If False, it
+        will not perform this check. Default is False.
 
     Examples
     ========
@@ -97,7 +102,9 @@ def DiscreteRV(symbol, density, set=S.Integers):
     set = sympify(set)
     pdf = Piecewise((density, set.as_relational(symbol)), (0, True))
     pdf = Lambda(symbol, pdf)
-    return rv(symbol.name, DiscreteDistributionHandmade, pdf, set)
+    # have a default of False while `rv` should have a default of True
+    kwargs['check'] = kwargs.pop('check', False)
+    return rv(symbol.name, DiscreteDistributionHandmade, pdf, set, **kwargs)
 
 
 #-------------------------------------------------------------------------------
diff --git a/sympy/stats/frv_types.py b/sympy/stats/frv_types.py
index ef89682832..5d2b7451ea 100644
--- a/sympy/stats/frv_types.py
+++ b/sympy/stats/frv_types.py
@@ -36,10 +36,11 @@
 'Rademacher'
 ]
 
-def rv(name, cls, *args):
+def rv(name, cls, *args, **kwargs):
     args = list(map(sympify, args))
     dist = cls(*args)
-    dist.check(*args)
+    if kwargs.pop('check', True):
+        dist.check(*args)
     pspace = SingleFinitePSpace(name, dist)
     if any(is_random(arg) for arg in args):
         from sympy.stats.compound_rv import CompoundPSpace, CompoundDistribution
@@ -66,17 +67,24 @@ def check(density):
         for p in density.values():
             _value_check((p >= 0, p <= 1),
                         "Probability at a point must be between 0 and 1.")
-        _value_check(Eq(sum(density.values()), 1), "Total Probability must be 1.")
+        val = sum(density.values())
+        _value_check(Eq(val, 1) != S.false, "Total Probability must be 1.")
 
-def FiniteRV(name, density):
+def FiniteRV(name, density, **kwargs):
     r"""
     Create a Finite Random Variable given a dict representing the density.
 
     Parameters
     ==========
 
+    name : Symbol
+        Represents name of the random variable.
     density: A dict
         Dictionary conatining the pdf of finite distribution
+    check : bool
+        If True, it will check whether the given density
+        integrates to 1 over the given set. If False, it
+        will not perform this check. Default is False.
 
     Examples
     ========
@@ -97,7 +105,9 @@ def FiniteRV(name, density):
     RandomSymbol
 
     """
-    return rv(name, FiniteDistributionHandmade, density)
+    # have a default of False while `rv` should have a default of True
+    kwargs['check'] = kwargs.pop('check', False)
+    return rv(name, FiniteDistributionHandmade, density, **kwargs)
 
 class DiscreteUniformDistribution(SingleFiniteDistribution):
 


Test patch:
diff --git a/sympy/stats/tests/test_continuous_rv.py b/sympy/stats/tests/test_continuous_rv.py
index 5b32b68878..16ce0d0e0f 100644
--- a/sympy/stats/tests/test_continuous_rv.py
+++ b/sympy/stats/tests/test_continuous_rv.py
@@ -337,7 +337,7 @@ def test_sample_continuous():
 def test_ContinuousRV():
     pdf = sqrt(2)*exp(-x**2/2)/(2*sqrt(pi))  # Normal distribution
     # X and Y should be equivalent
-    X = ContinuousRV(x, pdf)
+    X = ContinuousRV(x, pdf, check=True)
     Y = Normal('y', 0, 1)
 
     assert variance(X) == variance(Y)
@@ -346,7 +346,16 @@ def test_ContinuousRV():
     assert Z.pspace.domain.set == Interval(0, oo)
     assert E(Z) == 1
     assert P(Z > 5) == exp(-5)
-    raises(ValueError, lambda: ContinuousRV(z, exp(-z), set=Interval(0, 10)))
+    raises(ValueError, lambda: ContinuousRV(z, exp(-z), set=Interval(0, 10), check=True))
+
+    # the correct pdf for Gamma(k, theta) but the integral in `check`
+    # integrates to something equivalent to 1 and not to 1 exactly
+    _x, k, theta = symbols("x k theta", positive=True)
+    pdf = 1/(gamma(k)*theta**k)*_x**(k-1)*exp(-_x/theta)
+    X = ContinuousRV(_x, pdf, set=Interval(0, oo))
+    Y = Gamma('y', k, theta)
+    assert (E(X) - E(Y)).simplify() == 0
+    assert (variance(X) - variance(Y)).simplify() == 0
 
 
 def test_arcsin():
diff --git a/sympy/stats/tests/test_discrete_rv.py b/sympy/stats/tests/test_discrete_rv.py
index 30388483e6..4321587b42 100644
--- a/sympy/stats/tests/test_discrete_rv.py
+++ b/sympy/stats/tests/test_discrete_rv.py
@@ -178,11 +178,17 @@ def test_DiscreteRV():
     p = S(1)/2
     x = Symbol('x', integer=True, positive=True)
     pdf = p*(1 - p)**(x - 1) # pdf of Geometric Distribution
-    D = DiscreteRV(x, pdf, set=S.Naturals)
+    D = DiscreteRV(x, pdf, set=S.Naturals, check=True)
     assert E(D) == E(Geometric('G', S(1)/2)) == 2
     assert P(D > 3) == S(1)/8
     assert D.pspace.domain.set == S.Naturals
-    raises(ValueError, lambda: DiscreteRV(x, x, FiniteSet(*range(4))))
+    raises(ValueError, lambda: DiscreteRV(x, x, FiniteSet(*range(4)), check=True))
+
+    # purposeful invalid pmf but it should not raise since check=False
+    # see test_drv_types.test_ContinuousRV for explanation
+    X = DiscreteRV(x, 1/x, S.Naturals)
+    assert P(X < 2) == 1
+    assert E(X) == oo
 
 def test_precomputed_characteristic_functions():
     import mpmath
diff --git a/sympy/stats/tests/test_finite_rv.py b/sympy/stats/tests/test_finite_rv.py
index 5598df07b4..de2d825b41 100644
--- a/sympy/stats/tests/test_finite_rv.py
+++ b/sympy/stats/tests/test_finite_rv.py
@@ -383,7 +383,7 @@ def test_rademacher():
 
 
 def test_FiniteRV():
-    F = FiniteRV('F', {1: S.Half, 2: Rational(1, 4), 3: Rational(1, 4)})
+    F = FiniteRV('F', {1: S.Half, 2: Rational(1, 4), 3: Rational(1, 4)}, check=True)
     p = Symbol("p", positive=True)
 
     assert dict(density(F).items()) == {S.One: S.Half, S(2): Rational(1, 4), S(3): Rational(1, 4)}
@@ -395,10 +395,16 @@ def test_FiniteRV():
         *[Eq(F.symbol, i) for i in [1, 2, 3]])
 
     assert F.pspace.domain.set == FiniteSet(1, 2, 3)
-    raises(ValueError, lambda: FiniteRV('F', {1: S.Half, 2: S.Half, 3: S.Half}))
-    raises(ValueError, lambda: FiniteRV('F', {1: S.Half, 2: Rational(-1, 2), 3: S.One}))
+    raises(ValueError, lambda: FiniteRV('F', {1: S.Half, 2: S.Half, 3: S.Half}, check=True))
+    raises(ValueError, lambda: FiniteRV('F', {1: S.Half, 2: Rational(-1, 2), 3: S.One}, check=True))
     raises(ValueError, lambda: FiniteRV('F', {1: S.One, 2: Rational(3, 2), 3: S.Zero,\
-        4: Rational(-1, 2), 5: Rational(-3, 4), 6: Rational(-1, 4)}))
+        4: Rational(-1, 2), 5: Rational(-3, 4), 6: Rational(-1, 4)}, check=True))
+
+    # purposeful invalid pmf but it should not raise since check=False
+    # see test_drv_types.test_ContinuousRV for explanation
+    X = FiniteRV('X', {1: 1, 2: 2})
+    assert E(X) == 5
+    assert P(X <= 2) + P(X > 2) != 1
 
 def test_density_call():
     from sympy.abc import p