text: from sympy import E as e
  from sympy import (
      Symbol,
      Abs,
      exp,
      expint,
      S,
      pi,
      simplify,
      Interval,
      erf,
      erfc,
      Ne,
      EulerGamma,
      Eq,
      log,
      lowergamma,
      uppergamma,
      symbols,
      sqrt,
      And,
      gamma,
      beta,
      Piecewise,
      Integral,
      sin,
      cos,
      tan,
      sinh,
      cosh,
      besseli,
      floor,
      expand_func,
      Rational,
      I,
      re,
      im,
      lambdify,
      hyper,
      diff,
      Or,
      Mul,
      sign,
      Dummy,
      Sum,
      factorial,
      binomial,
      erfi,
      besselj,
  )
  from sympy.external import import_module
  from sympy.functions.special.error_functions import erfinv
  from sympy.functions.special.hyper import meijerg
  from sympy.sets.sets import Intersection, FiniteSet
  from sympy.stats import (
      P,
      E,
      where,
      density,
      variance,
      covariance,
      skewness,
      kurtosis,
      given,
      pspace,
      cdf,
      characteristic_function,
      moment_generating_function,
      ContinuousRV,
      sample,
      Arcsin,
      Benini,
      Beta,
      BetaNoncentral,
      BetaPrime,
      Cauchy,
      Chi,
      ChiSquared,
      ChiNoncentral,
      Dagum,
      Erlang,
      ExGaussian,
      Exponential,
      ExponentialPower,
      FDistribution,
      FisherZ,
      Frechet,
      Gamma,
      GammaInverse,
      Gompertz,
      Gumbel,
      Kumaraswamy,
      Laplace,
      Levy,
      Logistic,
      LogLogistic,
      LogNormal,
      Maxwell,
      Moyal,
      Nakagami,
      Normal,
      GaussianInverse,
      Pareto,
      PowerFunction,
      QuadraticU,
      RaisedCosine,
      Rayleigh,
      Reciprocal,
      ShiftedGompertz,
      StudentT,
      Trapezoidal,
      Triangular,
      Uniform,
      UniformSum,
      VonMises,
      Weibull,
      WignerSemicircle,
      Wald,
      correlation,
      moment,
      cmoment,
      smoment,
      quantile,
  )
  from sympy.stats.crv_types import NormalDistribution, ExponentialDistribution
  from sympy.stats.crv import SingleContinuousPSpace
  from sympy.stats.joint_rv import JointPSpace
  from sympy.testing.pytest import raises, XFAIL, slow, skip
  from sympy.testing.randtest import verify_numerically as tn
  
  oo = S.Infinity
  
  x, y, z = map(Symbol, "xyz")
  
  
  def test_single_normal():
      mu = Symbol("mu", real=True)
      sigma = Symbol("sigma", positive=True)
      X = Normal("x", 0, 1)
      Y = X * sigma + mu
  
      assert E(Y) == mu
      assert variance(Y) == sigma ** 2
      pdf = density(Y)
      x = Symbol("x", real=True)
      assert pdf(x) == 2 ** S.Half * exp(-((x - mu) ** 2) / (2 * sigma ** 2)) / (
          2 * pi ** S.Half * sigma
      )
  
      assert P(X ** 2 < 1) == erf(2 ** S.Half / 2)
      assert quantile(Y)(x) == Intersection(
          S.Reals,
          FiniteSet(sqrt(2) * sigma * (sqrt(2) * mu / (2 * sigma) + erfinv(2 * x - 1))),
      )
      assert E(X, Eq(X, mu)) == mu
  
      # issue 8248
      assert X.pspace.compute_expectation(1).doit() == 1
  
  
  def test_conditional_1d():
      X = Normal("x", 0, 1)
      Y = given(X, X >= 0)
      z = Symbol("z")
  
      assert density(Y)(z) == 2 * density(X)(z)
  
      assert Y.pspace.domain.set == Interval(0, oo)
      assert E(Y) == sqrt(2) / sqrt(pi)
  
      assert E(X ** 2) == E(Y ** 2)
  
  
  def test_ContinuousDomain():
      X = Normal("x", 0, 1)
      assert where(X ** 2 <= 1).set == Interval(-1, 1)
      assert where(X ** 2 <= 1).symbol == X.symbol
      where(And(X ** 2 <= 1, X >= 0)).set == Interval(0, 1)
      raises(ValueError, lambda: where(sin(X) > 1))
  
      Y = given(X, X >= 0)
  
      assert Y.pspace.domain.set == Interval(0, oo)
  
  
  @slow
  def test_multiple_normal():
      X, Y = Normal("x", 0, 1), Normal("y", 0, 1)
      p = Symbol("p", positive=True)
  
      assert E(X + Y) == 0
      assert variance(X + Y) == 2
      assert variance(X + X) == 4
      assert covariance(X, Y) == 0
      assert covariance(2 * X + Y, -X) == -2 * variance(X)
      assert skewness(X) == 0
      assert skewness(X + Y) == 0
      assert kurtosis(X) == 3
      assert kurtosis(X + Y) == 3
      assert correlation(X, Y) == 0
      assert correlation(X, X + Y) == correlation(X, X - Y)
      assert moment(X, 2) == 1
      assert cmoment(X, 3) == 0
      assert moment(X + Y, 4) == 12
      assert cmoment(X, 2) == variance(X)
      assert smoment(X * X, 2) == 1
      assert smoment(X + Y, 3) == skewness(X + Y)
      assert smoment(X + Y, 4) == kurtosis(X + Y)
      assert E(X, Eq(X + Y, 0)) == 0
      assert variance(X, Eq(X + Y, 0)) == S.Half
      assert quantile(X)(p) == sqrt(2) * erfinv(2 * p - S.One)
  
  
  def test_symbolic():
      mu1, mu2 = symbols("mu1 mu2", real=True)
      s1, s2 = symbols("sigma1 sigma2", positive=True)
      rate = Symbol("lambda", positive=True)
      X = Normal("x", mu1, s1)
      Y = Normal("y", mu2, s2)
      Z = Exponential("z", rate)
      a, b, c = symbols("a b c", real=True)
  
      assert E(X) == mu1
      assert E(X + Y) == mu1 + mu2
      assert E(a * X + b) == a * E(X) + b
      assert variance(X) == s1 ** 2
      assert variance(X + a * Y + b) == variance(X) + a ** 2 * variance(Y)
  
      assert E(Z) == 1 / rate
      assert E(a * Z + b) == a * E(Z) + b
      assert E(X + a * Z + b) == mu1 + a / rate + b
  
  
  def test_cdf():
      X = Normal("x", 0, 1)
  
      d = cdf(X)
      assert P(X < 1) == d(1).rewrite(erfc)
      assert d(0) == S.Half
  
      d = cdf(X, X > 0)  # given X>0
      assert d(0) == 0
  
      Y = Exponential("y", 10)
      d = cdf(Y)
      assert d(-5) == 0
      assert P(Y > 3) == 1 - d(3)
  
      raises(ValueError, lambda: cdf(X + Y))
  
      Z = Exponential("z", 1)
      f = cdf(Z)
      assert f(z) == Piecewise((1 - exp(-z), z >= 0), (0, True))
  
  
  def test_characteristic_function():
      X = Uniform("x", 0, 1)
  
      cf = characteristic_function(X)
      assert cf(1) == -I * (-1 + exp(I))
  
      Y = Normal("y", 1, 1)
      cf = characteristic_function(Y)
      assert cf(0) == 1
      assert cf(1) == exp(I - S.Half)
  
      Z = Exponential("z", 5)
      cf = characteristic_function(Z)
      assert cf(0) == 1
      assert cf(1).expand() == Rational(25, 26) + I * Rational(5, 26)
  
      X = GaussianInverse("x", 1, 1)
      cf = characteristic_function(X)
      assert cf(0) == 1
      assert cf(1) == exp(1 - sqrt(1 - 2 * I))
  
      X = ExGaussian("x", 0, 1, 1)
      cf = characteristic_function(X)
      assert cf(0) == 1
      assert cf(1) == (1 + I) * exp(Rational(-1, 2)) / 2
  
      L = Levy("x", 0, 1)
      cf = characteristic_function(L)
      assert cf(0) == 1
      assert cf(1) == exp(-sqrt(2) * sqrt(-I))
  
  
  def test_moment_generating_function():
      t = symbols("t", positive=True)
  
      # Symbolic tests
      a, b, c = symbols("a b c")
  
      mgf = moment_generating_function(Beta("x", a, b))(t)
      assert mgf == hyper((a,), (a + b,), t)
  
      mgf = moment_generating_function(Chi("x", a))(t)
      assert mgf == sqrt(2) * t * gamma(a / 2 + S.Half) * hyper(
          (a / 2 + S.Half,), (Rational(3, 2),), t ** 2 / 2
      ) / gamma(a / 2) + hyper((a / 2,), (S.Half,), t ** 2 / 2)
  
      mgf = moment_generating_function(ChiSquared("x", a))(t)
      assert mgf == (1 - 2 * t) ** (-a / 2)
  
      mgf = moment_generating_function(Erlang("x", a, b))(t)
      assert mgf == (1 - t / b) ** (-a)
  
      mgf = moment_generating_function(ExGaussian("x", a, b, c))(t)
      assert mgf == exp(a * t + b ** 2 * t ** 2 / 2) / (1 - t / c)
  
      mgf = moment_generating_function(Exponential("x", a))(t)
      assert mgf == a / (a - t)
  
      mgf = moment_generating_function(Gamma("x", a, b))(t)
      assert mgf == (-b * t + 1) ** (-a)
  
      mgf = moment_generating_function(Gumbel("x", a, b))(t)
      assert mgf == exp(b * t) * gamma(-a * t + 1)
  
      mgf = moment_generating_function(Gompertz("x", a, b))(t)
      assert mgf == b * exp(b) * expint(t / a, b)
  
      mgf = moment_generating_function(Laplace("x", a, b))(t)
      assert mgf == exp(a * t) / (-(b ** 2) * t ** 2 + 1)
  
      mgf = moment_generating_function(Logistic("x", a, b))(t)
      assert mgf == exp(a * t) * beta(-b * t + 1, b * t + 1)
  
      mgf = moment_generating_function(Normal("x", a, b))(t)
      assert mgf == exp(a * t + b ** 2 * t ** 2 / 2)
  
      mgf = moment_generating_function(Pareto("x", a, b))(t)
      assert mgf == b * (-a * t) ** b * uppergamma(-b, -a * t)
  
      mgf = moment_generating_function(QuadraticU("x", a, b))(t)
      assert str(mgf) == (
          "(3*(t*(-4*b + (a + b)**2) + 4)*exp(b*t) - "
          "3*(t*(a**2 + 2*a*(b - 2) + b**2) + 4)*exp(a*t))/(t**2*(a - b)**3)"
      )
  
      mgf = moment_generating_function(RaisedCosine("x", a, b))(t)
      assert mgf == pi ** 2 * exp(a * t) * sinh(b * t) / (
          b * t * (b ** 2 * t ** 2 + pi ** 2)
      )
  
      mgf = moment_generating_function(Rayleigh("x", a))(t)
      assert (
          mgf
          == sqrt(2)
          * sqrt(pi)
          * a
          * t
          * (erf(sqrt(2) * a * t / 2) + 1)
          * exp(a ** 2 * t ** 2 / 2)
          / 2
          + 1
      )
  
      mgf = moment_generating_function(Triangular("x", a, b, c))(t)
      assert str(mgf) == (
          "(-2*(-a + b)*exp(c*t) + 2*(-a + c)*exp(b*t) + "
          "2*(b - c)*exp(a*t))/(t**2*(-a + b)*(-a + c)*(b - c))"
      )
  
      mgf = moment_generating_function(Uniform("x", a, b))(t)
      assert mgf == (-exp(a * t) + exp(b * t)) / (t * (-a + b))
  
      mgf = moment_generating_function(UniformSum("x", a))(t)
      assert mgf == ((exp(t) - 1) / t) ** a
  
      mgf = moment_generating_function(WignerSemicircle("x", a))(t)
      assert mgf == 2 * besseli(1, a * t) / (a * t)
  
      # Numeric tests
  
      mgf = moment_generating_function(Beta("x", 1, 1))(t)
      assert mgf.diff(t).subs(t, 1) == hyper((2,), (3,), 1) / 2
  
      mgf = moment_generating_function(Chi("x", 1))(t)
      assert mgf.diff(t).subs(t, 1) == sqrt(2) * hyper(
          (1,), (Rational(3, 2),), S.Half
      ) / sqrt(pi) + hyper((Rational(3, 2),), (Rational(3, 2),), S.Half) + 2 * sqrt(
          2
      ) * hyper(
          (2,), (Rational(5, 2),), S.Half
      ) / (
          3 * sqrt(pi)
      )
  
      mgf = moment_generating_function(ChiSquared("x", 1))(t)
      assert mgf.diff(t).subs(t, 1) == I
  
      mgf = moment_generating_function(Erlang("x", 1, 1))(t)
      assert mgf.diff(t).subs(t, 0) == 1
  
      mgf = moment_generating_function(ExGaussian("x", 0, 1, 1))(t)
      assert mgf.diff(t).subs(t, 2) == -exp(2)
  
      mgf = moment_generating_function(Exponential("x", 1))(t)
      assert mgf.diff(t).subs(t, 0) == 1
  
      mgf = moment_generating_function(Gamma("x", 1, 1))(t)
      assert mgf.diff(t).subs(t, 0) == 1
  
      mgf = moment_generating_function(Gumbel("x", 1, 1))(t)
      assert mgf.diff(t).subs(t, 0) == EulerGamma + 1
  
      mgf = moment_generating_function(Gompertz("x", 1, 1))(t)
      assert mgf.diff(t).subs(t, 1) == -e * meijerg(((), (1, 1)), ((0, 0, 0), ()), 1)
  
      mgf = moment_generating_function(Laplace("x", 1, 1))(t)
      assert mgf.diff(t).subs(t, 0) == 1
  
      mgf = moment_generating_function(Logistic("x", 1, 1))(t)
      assert mgf.diff(t).subs(t, 0) == beta(1, 1)
  
      mgf = moment_generating_function(Normal("x", 0, 1))(t)
      assert mgf.diff(t).subs(t, 1) == exp(S.Half)
  
      mgf = moment_generating_function(Pareto("x", 1, 1))(t)
      assert mgf.diff(t).subs(t, 0) == expint(1, 0)
  
      mgf = moment_generating_function(QuadraticU("x", 1, 2))(t)
      assert mgf.diff(t).subs(t, 1) == -12 * e - 3 * exp(2)
  
      mgf = moment_generating_function(RaisedCosine("x", 1, 1))(t)
      assert mgf.diff(t).subs(t, 1) == -2 * e * pi ** 2 * sinh(1) / (
          1 + pi ** 2
      ) ** 2 + e * pi ** 2 * cosh(1) / (1 + pi ** 2)
  
      mgf = moment_generating_function(Rayleigh("x", 1))(t)
      assert mgf.diff(t).subs(t, 0) == sqrt(2) * sqrt(pi) / 2
  
      mgf = moment_generating_function(Triangular("x", 1, 3, 2))(t)
      assert mgf.diff(t).subs(t, 1) == -e + exp(3)
  
      mgf = moment_generating_function(Uniform("x", 0, 1))(t)
      assert mgf.diff(t).subs(t, 1) == 1
  
      mgf = moment_generating_function(UniformSum("x", 1))(t)
      assert mgf.diff(t).subs(t, 1) == 1
  
      mgf = moment_generating_function(WignerSemicircle("x", 1))(t)
      assert mgf.diff(t).subs(t, 1) == -2 * besseli(1, 1) + besseli(2, 1) + besseli(0, 1)
  
  
  def test_sample_continuous():
      Z = ContinuousRV(z, exp(-z), set=Interval(0, oo))
      assert sample(Z) in Z.pspace.domain.set
      sym, val = list(Z.pspace.sample().items())[0]
      assert sym == Z and val in Interval(0, oo)
      assert density(Z)(-1) == 0
  
  
  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)
      Y = Normal("y", 0, 1)
  
      assert variance(X) == variance(Y)
      assert P(X > 0) == P(Y > 0)
  
  
  def test_arcsin():
      from sympy import asin
  
      a = Symbol("a", real=True)
      b = Symbol("b", real=True)
  
      X = Arcsin("x", a, b)
      assert density(X)(x) == 1 / (pi * sqrt((-x + b) * (x - a)))
      assert cdf(X)(x) == Piecewise(
          (0, a > x), (2 * asin(sqrt((-a + x) / (-a + b))) / pi, b >= x), (1, True)
      )
  
  
  def test_benini():
      alpha = Symbol("alpha", positive=True)
      beta = Symbol("beta", positive=True)
      sigma = Symbol("sigma", positive=True)
      X = Benini("x", alpha, beta, sigma)
  
      assert density(X)(x) == (
          (alpha / x + 2 * beta * log(x / sigma) / x)
          * exp(-alpha * log(x / sigma) - beta * log(x / sigma) ** 2)
      )
  
      alpha = Symbol("alpha", nonpositive=True)
      raises(ValueError, lambda: Benini("x", alpha, beta, sigma))
  
      beta = Symbol("beta", nonpositive=True)
      raises(ValueError, lambda: Benini("x", alpha, beta, sigma))
  
      alpha = Symbol("alpha", positive=True)
      raises(ValueError, lambda: Benini("x", alpha, beta, sigma))
  
      beta = Symbol("beta", positive=True)
      sigma = Symbol("sigma", nonpositive=True)
      raises(ValueError, lambda: Benini("x", alpha, beta, sigma))
  
  
  def test_beta():
      a, b = symbols("alpha beta", positive=True)
      B = Beta("x", a, b)
  
      assert pspace(B).domain.set == Interval(0, 1)
      assert characteristic_function(B)(x) == hyper((a,), (a + b,), I * x)
      assert density(B)(x) == x ** (a - 1) * (1 - x) ** (b - 1) / beta(a, b)
  
      assert simplify(E(B)) == a / (a + b)
      assert simplify(variance(B)) == a * b / (
          a ** 3 + 3 * a ** 2 * b + a ** 2 + 3 * a * b ** 2 + 2 * a * b + b ** 3 + b ** 2
      )
  
      # Full symbolic solution is too much, test with numeric version
      a, b = 1, 2
      B = Beta("x", a, b)
      assert expand_func(E(B)) == a / S(a + b)
      assert expand_func(variance(B)) == (a * b) / S((a + b) ** 2 * (a + b + 1))
  
  
  def test_beta_noncentral():
      a, b = symbols("a b", positive=True)
      c = Symbol("c", nonnegative=True)
      _k = Dummy("k")
  
      X = BetaNoncentral("x", a, b, c)
  
      assert pspace(X).domain.set == Interval(0, 1)
  
      dens = density(X)
      z = Symbol("z")
  
      res = Sum(
          z ** (_k + a - 1)
          * (c / 2) ** _k
          * (1 - z) ** (b - 1)
          * exp(-c / 2)
          / (beta(_k + a, b) * factorial(_k)),
          (_k, 0, oo),
      )
      assert dens(z).dummy_eq(res)
  
      # BetaCentral should not raise if the assumptions
      # on the symbols can not be determined
      a, b, c = symbols("a b c")
      assert BetaNoncentral("x", a, b, c)
  
      a = Symbol("a", positive=False, real=True)
      raises(ValueError, lambda: BetaNoncentral("x", a, b, c))
  
      a = Symbol("a", positive=True)
      b = Symbol("b", positive=False, real=True)
      raises(ValueError, lambda: BetaNoncentral("x", a, b, c))
  
      a = Symbol("a", positive=True)
      b = Symbol("b", positive=True)
      c = Symbol("c", nonnegative=False, real=True)
      raises(ValueError, lambda: BetaNoncentral("x", a, b, c))
  
  
  def test_betaprime():
      alpha = Symbol("alpha", positive=True)
  
      betap = Symbol("beta", positive=True)
  
      X = BetaPrime("x", alpha, betap)
      assert density(X)(x) == x ** (alpha - 1) * (x + 1) ** (-alpha - betap) / beta(
          alpha, betap
      )
  
      alpha = Symbol("alpha", nonpositive=True)
      raises(ValueError, lambda: BetaPrime("x", alpha, betap))
  
      alpha = Symbol("alpha", positive=True)
      betap = Symbol("beta", nonpositive=True)
      raises(ValueError, lambda: BetaPrime("x", alpha, betap))
  
  
  def test_cauchy():
      x0 = Symbol("x0")
      gamma = Symbol("gamma", positive=True)
      p = Symbol("p", positive=True)
  
      X = Cauchy("x", x0, gamma)
      # Tests the characteristic function
      assert characteristic_function(X)(x) == exp(-gamma * Abs(x) + I * x * x0)
  
      assert density(X)(x) == 1 / (pi * gamma * (1 + (x - x0) ** 2 / gamma ** 2))
      assert diff(cdf(X)(x), x) == density(X)(x)
      assert quantile(X)(p) == gamma * tan(pi * (p - S.Half)) + x0
  
      gamma = Symbol("gamma", nonpositive=True)
      raises(ValueError, lambda: Cauchy("x", x0, gamma))
  
  
  def test_chi():
      from sympy import I
  
      k = Symbol("k", integer=True)
  
      X = Chi("x", k)
      assert density(X)(x) == 2 ** (-k / 2 + 1) * x ** (k - 1) * exp(
          -(x ** 2) / 2
      ) / gamma(k / 2)
  
      # Tests the characteristic function
      assert characteristic_function(X)(x) == sqrt(2) * I * x * gamma(
          k / 2 + S(1) / 2
      ) * hyper((k / 2 + S(1) / 2,), (S(3) / 2,), -(x ** 2) / 2) / gamma(k / 2) + hyper(
          (k / 2,), (S(1) / 2,), -(x ** 2) / 2
      )
  
      # Tests the moment generating function
      assert moment_generating_function(X)(x) == sqrt(2) * x * gamma(
          k / 2 + S(1) / 2
      ) * hyper((k / 2 + S(1) / 2,), (S(3) / 2,), x ** 2 / 2) / gamma(k / 2) + hyper(
          (k / 2,), (S(1) / 2,), x ** 2 / 2
      )
  
      k = Symbol("k", integer=True, positive=False)
      raises(ValueError, lambda: Chi("x", k))
  
      k = Symbol("k", integer=False, positive=True)
      raises(ValueError, lambda: Chi("x", k))
  
  
  def test_chi_noncentral():
      k = Symbol("k", integer=True)
      l = Symbol("l")
  
      X = ChiNoncentral("x", k, l)
      assert density(X)(x) == (
          x ** k
          * l
          * (x * l) ** (-k / 2)
          * exp(-(x ** 2) / 2 - l ** 2 / 2)
          * besseli(k / 2 - 1, x * l)
      )
  
      k = Symbol("k", integer=True, positive=False)
      raises(ValueError, lambda: ChiNoncentral("x", k, l))
  
      k = Symbol("k", integer=True, positive=True)
      l = Symbol("l", nonpositive=True)
      raises(ValueError, lambda: ChiNoncentral("x", k, l))
  
      k = Symbol("k", integer=False)
      l = Symbol("l", positive=True)
      raises(ValueError, lambda: ChiNoncentral("x", k, l))
  
  
  def test_chi_squared():
      k = Symbol("k", integer=True)
      X = ChiSquared("x", k)
  
      # Tests the characteristic function
      assert characteristic_function(X)(x) == ((-2 * I * x + 1) ** (-k / 2))
  
      assert density(X)(x) == 2 ** (-k / 2) * x ** (k / 2 - 1) * exp(-x / 2) / gamma(
          k / 2
      )
      assert cdf(X)(x) == Piecewise(
          (lowergamma(k / 2, x / 2) / gamma(k / 2), x >= 0), (0, True)
      )
      assert E(X) == k
      assert variance(X) == 2 * k
  
      X = ChiSquared("x", 15)
      assert cdf(X)(3) == -14873 * sqrt(6) * exp(Rational(-3, 2)) / (
          5005 * sqrt(pi)
      ) + erf(sqrt(6) / 2)
  
      k = Symbol("k", integer=True, positive=False)
      raises(ValueError, lambda: ChiSquared("x", k))
  
      k = Symbol("k", integer=False, positive=True)
      raises(ValueError, lambda: ChiSquared("x", k))
  
  
  def test_dagum():
      p = Symbol("p", positive=True)
      b = Symbol("b", positive=True)
      a = Symbol("a", positive=True)
  
      X = Dagum("x", p, a, b)
      assert (
          density(X)(x) == a * p * (x / b) ** (a * p) * ((x / b) ** a + 1) ** (-p - 1) / x
      )
      assert cdf(X)(x) == Piecewise(((1 + (x / b) ** (-a)) ** (-p), x >= 0), (0, True))
  
      p = Symbol("p", nonpositive=True)
      raises(ValueError, lambda: Dagum("x", p, a, b))
  
      p = Symbol("p", positive=True)
      b = Symbol("b", nonpositive=True)
      raises(ValueError, lambda: Dagum("x", p, a, b))
  
      b = Symbol("b", positive=True)
      a = Symbol("a", nonpositive=True)
      raises(ValueError, lambda: Dagum("x", p, a, b))
  
  
  def test_erlang():
      k = Symbol("k", integer=True, positive=True)
      l = Symbol("l", positive=True)
  
      X = Erlang("x", k, l)
      assert density(X)(x) == x ** (k - 1) * l ** k * exp(-x * l) / gamma(k)
      assert cdf(X)(x) == Piecewise((lowergamma(k, l * x) / gamma(k), x > 0), (0, True))
  
  
  def test_exgaussian():
      m, z = symbols("m, z")
      s, l = symbols("s, l", positive=True)
      X = ExGaussian("x", m, s, l)
  
      assert (
          density(X)(z)
          == l
          * exp(l * (l * s ** 2 + 2 * m - 2 * z) / 2)
          * erfc(sqrt(2) * (l * s ** 2 + m - z) / (2 * s))
          / 2
      )
  
      # Note: actual_output simplifies to expected_output.
      # Ideally cdf(X)(z) would return expected_output
      # expected_output = (erf(sqrt(2)*(l*s**2 + m - z)/(2*s)) - 1)*exp(l*(l*s**2 + 2*m - 2*z)/2)/2 - erf(sqrt(2)*(m - z)/(2*s))/2 + S.Half
      u = l * (z - m)
      v = l * s
      GaussianCDF1 = cdf(Normal("x", 0, v))(u)
      GaussianCDF2 = cdf(Normal("x", v ** 2, v))(u)
      actual_output = GaussianCDF1 - exp(-u + (v ** 2 / 2) + log(GaussianCDF2))
      assert cdf(X)(z) == actual_output
      # assert simplify(actual_output) == expected_output
  
      assert variance(X).expand() == s ** 2 + l ** (-2)
  
      assert skewness(X).expand() == 2 / (
          l ** 3 * s ** 2 * sqrt(s ** 2 + l ** (-2)) + l * sqrt(s ** 2 + l ** (-2))
      )
  
  
  def test_exponential():
      rate = Symbol("lambda", positive=True)
      X = Exponential("x", rate)
      p = Symbol("p", positive=True, real=True, finite=True)
  
      assert E(X) == 1 / rate
      assert variance(X) == 1 / rate ** 2
      assert skewness(X) == 2
      assert skewness(X) == smoment(X, 3)
      assert kurtosis(X) == 9
      assert kurtosis(X) == smoment(X, 4)
      assert smoment(2 * X, 4) == smoment(X, 4)
      assert moment(X, 3) == 3 * 2 * 1 / rate ** 3
      assert P(X > 0) is S.One
      assert P(X > 1) == exp(-rate)
      assert P(X > 10) == exp(-10 * rate)
      assert quantile(X)(p) == -log(1 - p) / rate
  
      assert where(X <= 1).set == Interval(0, 1)
      # Test issue 9970
      z = Dummy("z")
      assert P(X > z) == exp(-z * rate)
      assert P(X < z) == 0
      # Test issue 10076 (Distribution with interval(0,oo))
      x = Symbol("x")
      _z = Dummy("_z")
      b = SingleContinuousPSpace(x, ExponentialDistribution(2))
  
      expected1 = Integral(2 * exp(-2 * _z), (_z, 3, oo))
      assert b.probability(x > 3, evaluate=False).dummy_eq(expected1) is True
  
      expected2 = Integral(2 * exp(-2 * _z), (_z, 0, 4))
      assert b.probability(x < 4, evaluate=False).dummy_eq(expected2) is True
  
  
  def test_exponential_power():
      mu = Symbol("mu")
      z = Symbol("z")
      alpha = Symbol("alpha", positive=True)
      beta = Symbol("beta", positive=True)
  
      X = ExponentialPower("x", mu, alpha, beta)
  
      assert density(X)(z) == beta * exp(-((Abs(mu - z) / alpha) ** beta)) / (
          2 * alpha * gamma(1 / beta)
      )
      assert cdf(X)(z) == S.Half + lowergamma(
          1 / beta, (Abs(mu - z) / alpha) ** beta
      ) * sign(-mu + z) / (2 * gamma(1 / beta))
  
  
  def test_f_distribution():
      d1 = Symbol("d1", positive=True)
      d2 = Symbol("d2", positive=True)
  
      X = FDistribution("x", d1, d2)
  
      assert density(X)(x) == (
          d2 ** (d2 / 2)
          * sqrt((d1 * x) ** d1 * (d1 * x + d2) ** (-d1 - d2))
          / (x * beta(d1 / 2, d2 / 2))
      )
  
      d1 = Symbol("d1", nonpositive=True)
      raises(ValueError, lambda: FDistribution("x", d1, d1))
  
      d1 = Symbol("d1", positive=True, integer=False)
      raises(ValueError, lambda: FDistribution("x", d1, d1))
  
      d1 = Symbol("d1", positive=True)
      d2 = Symbol("d2", nonpositive=True)
      raises(ValueError, lambda: FDistribution("x", d1, d2))
  
      d2 = Symbol("d2", positive=True, integer=False)
      raises(ValueError, lambda: FDistribution("x", d1, d2))
  
  
  def test_fisher_z():
      d1 = Symbol("d1", positive=True)
      d2 = Symbol("d2", positive=True)
  
      X = FisherZ("x", d1, d2)
      assert density(X)(x) == (
          2
          * d1 ** (d1 / 2)
          * d2 ** (d2 / 2)
          * (d1 * exp(2 * x) + d2) ** (-d1 / 2 - d2 / 2)
          * exp(d1 * x)
          / beta(d1 / 2, d2 / 2)
      )
  
  
  def test_frechet():
      a = Symbol("a", positive=True)
      s = Symbol("s", positive=True)
      m = Symbol("m", real=True)
  
      X = Frechet("x", a, s=s, m=m)
      assert (
          density(X)(x)
          == a * ((x - m) / s) ** (-a - 1) * exp(-(((x - m) / s) ** (-a))) / s
      )
      assert cdf(X)(x) == Piecewise((exp(-(((-m + x) / s) ** (-a))), m <= x), (0, True))
  
  
  def test_gamma():
      k = Symbol("k", positive=True)
      theta = Symbol("theta", positive=True)
  
      X = Gamma("x", k, theta)
  
      # Tests characteristic function
      assert characteristic_function(X)(x) == ((-I * theta * x + 1) ** (-k))
  
      assert density(X)(x) == x ** (k - 1) * theta ** (-k) * exp(-x / theta) / gamma(k)
      assert cdf(X, meijerg=True)(z) == Piecewise(
          (
              -k * lowergamma(k, 0) / gamma(k + 1)
              + k * lowergamma(k, z / theta) / gamma(k + 1),
              z >= 0,
          ),
          (0, True),
      )
  
      # assert simplify(variance(X)) == k*theta**2  # handled numerically below
      assert E(X) == moment(X, 1)
  
      k, theta = symbols("k theta", positive=True)
      X = Gamma("x", k, theta)
      assert E(X) == k * theta
      assert variance(X) == k * theta ** 2
      assert skewness(X).expand() == 2 / sqrt(k)
      assert kurtosis(X).expand() == 3 + 6 / k
  
  
  def test_gamma_inverse():
      a = Symbol("a", positive=True)
      b = Symbol("b", positive=True)
      X = GammaInverse("x", a, b)
      assert density(X)(x) == x ** (-a - 1) * b ** a * exp(-b / x) / gamma(a)
      assert cdf(X)(x) == Piecewise((uppergamma(a, b / x) / gamma(a), x > 0), (0, True))
  
  
  def test_sampling_gamma_inverse():
      scipy = import_module("scipy")
      if not scipy:
          skip("Scipy not installed. Abort tests for sampling of gamma inverse.")
      X = GammaInverse("x", 1, 1)
      assert sample(X) in X.pspace.domain.set
  
  
  def test_gompertz():
      b = Symbol("b", positive=True)
      eta = Symbol("eta", positive=True)
  
      X = Gompertz("x", b, eta)
  
      assert density(X)(x) == b * eta * exp(eta) * exp(b * x) * exp(-eta * exp(b * x))
      assert cdf(X)(x) == 1 - exp(eta) * exp(-eta * exp(b * x))
      assert diff(cdf(X)(x), x) == density(X)(x)
  
  
  def test_gumbel():
      beta = Symbol("beta", positive=True)
      mu = Symbol("mu")
      x = Symbol("x")
      y = Symbol("y")
      X = Gumbel("x", beta, mu)
      Y = Gumbel("y", beta, mu, minimum=True)
      assert (
          density(X)(x).expand()
          == exp(mu / beta)
          * exp(-x / beta)
          * exp(-exp(mu / beta) * exp(-x / beta))
          / beta
      )
      assert (
          density(Y)(y).expand()
          == exp(-mu / beta)
          * exp(y / beta)
          * exp(-exp(-mu / beta) * exp(y / beta))
          / beta
      )
      assert cdf(X)(x).expand() == exp(-exp(mu / beta) * exp(-x / beta))
  
  
  def test_kumaraswamy():
      a = Symbol("a", positive=True)
      b = Symbol("b", positive=True)
  
      X = Kumaraswamy("x", a, b)
      assert density(X)(x) == x ** (a - 1) * a * b * (-(x ** a) + 1) ** (b - 1)
      assert cdf(X)(x) == Piecewise(
          (0, x < 0), (-((-(x ** a) + 1) ** b) + 1, x <= 1), (1, True)
      )
  
  
  def test_laplace():
      mu = Symbol("mu")
      b = Symbol("b", positive=True)
  
      X = Laplace("x", mu, b)
  
      # Tests characteristic_function
      assert characteristic_function(X)(x) == (exp(I * mu * x) / (b ** 2 * x ** 2 + 1))
  
      assert density(X)(x) == exp(-Abs(x - mu) / b) / (2 * b)
      assert cdf(X)(x) == Piecewise(
          (exp((-mu + x) / b) / 2, mu > x), (-exp((mu - x) / b) / 2 + 1, True)
      )
  
  
  def test_levy():
      mu = Symbol("mu", real=True)
      c = Symbol("c", positive=True)
  
      X = Levy("x", mu, c)
      assert X.pspace.domain.set == Interval(mu, oo)
      assert density(X)(x) == sqrt(c / (2 * pi)) * exp(-c / (2 * (x - mu))) / (
          (x - mu) ** (S.One + S.Half)
      )
      assert cdf(X)(x) == erfc(sqrt(c / (2 * (x - mu))))
  
      mu = Symbol("mu", real=False)
      raises(ValueError, lambda: Levy("x", mu, c))
  
      c = Symbol("c", nonpositive=True)
      raises(ValueError, lambda: Levy("x", mu, c))
  
      mu = Symbol("mu", real=True)
      raises(ValueError, lambda: Levy("x", mu, c))
  
  
  def test_logistic():
      mu = Symbol("mu", real=True)
      s = Symbol("s", positive=True)
      p = Symbol("p", positive=True)
  
      X = Logistic("x", mu, s)
  
      # Tests characteristics_function
      assert characteristic_function(X)(x) == (
          Piecewise(
              (pi * s * x * exp(I * mu * x) / sinh(pi * s * x), Ne(x, 0)), (1, True)
          )
      )
  
      assert density(X)(x) == exp((-x + mu) / s) / (s * (exp((-x + mu) / s) + 1) ** 2)
      assert cdf(X)(x) == 1 / (exp((mu - x) / s) + 1)
      assert quantile(X)(p) == mu - s * log(-S.One + 1 / p)
  
  
  def test_loglogistic():
      a, b = symbols("a b")
      assert LogLogistic("x", a, b)
  
      a = Symbol("a", negative=True)
      b = Symbol("b", positive=True)
      raises(ValueError, lambda: LogLogistic("x", a, b))
  
      a = Symbol("a", positive=True)
      b = Symbol("b", negative=True)
      raises(ValueError, lambda: LogLogistic("x", a, b))
  
      a, b, z, p = symbols("a b z p", positive=True)
      X = LogLogistic("x", a, b)
      assert density(X)(z) == b * (z / a) ** (b - 1) / (a * ((z / a) ** b + 1) ** 2)
      assert cdf(X)(z) == 1 / (1 + (z / a) ** (-b))
      assert quantile(X)(p) == a * (p / (1 - p)) ** (1 / b)
  
      # Expectation
      assert E(X) == Piecewise((S.NaN, b <= 1), (pi * a / (b * sin(pi / b)), True))
      b = symbols("b", prime=True)  # b > 1
      X = LogLogistic("x", a, b)
      assert E(X) == pi * a / (b * sin(pi / b))
  
  
  def test_lognormal():
      mean = Symbol("mu", real=True)
      std = Symbol("sigma", positive=True)
      X = LogNormal("x", mean, std)
      # The sympy integrator can't do this too well
      # assert E(X) == exp(mean+std**2/2)
      # assert variance(X) == (exp(std**2)-1) * exp(2*mean + std**2)
  
      # Right now, only density function and sampling works
  
      for i in range(3):
          X = LogNormal("x", i, 1)
          assert sample(X) in X.pspace.domain.set
      # The sympy integrator can't do this too well
      # assert E(X) ==
  
      mu = Symbol("mu", real=True)
      sigma = Symbol("sigma", positive=True)
  
      X = LogNormal("x", mu, sigma)
      assert density(X)(x) == (
          sqrt(2)
          * exp(-((-mu + log(x)) ** 2) / (2 * sigma ** 2))
          / (2 * x * sqrt(pi) * sigma)
      )
      # Tests cdf
      assert cdf(X)(x) == Piecewise(
          (erf(sqrt(2) * (-mu + log(x)) / (2 * sigma)) / 2 + S(1) / 2, x > 0), (0, True)
      )
  
      X = LogNormal("x", 0, 1)  # Mean 0, standard deviation 1
      assert density(X)(x) == sqrt(2) * exp(-log(x) ** 2 / 2) / (2 * x * sqrt(pi))
  
  
  def test_maxwell():
      a = Symbol("a", positive=True)
  
      X = Maxwell("x", a)
  
      assert density(X)(x) == (
          sqrt(2) * x ** 2 * exp(-(x ** 2) / (2 * a ** 2)) / (sqrt(pi) * a ** 3)
      )
      assert E(X) == 2 * sqrt(2) * a / sqrt(pi)
      assert variance(X) == -8 * a ** 2 / pi + 3 * a ** 2
      assert cdf(X)(x) == erf(sqrt(2) * x / (2 * a)) - sqrt(2) * x * exp(
          -(x ** 2) / (2 * a ** 2)
      ) / (sqrt(pi) * a)
      assert diff(cdf(X)(x), x) == density(X)(x)
  
  
  def test_Moyal():
      mu = Symbol("mu", real=False)
      sigma = Symbol("sigma", real=True, positive=True)
      raises(ValueError, lambda: Moyal("M", mu, sigma))
  
      mu = Symbol("mu", real=True)
      sigma = Symbol("sigma", real=True, negative=True)
      raises(ValueError, lambda: Moyal("M", mu, sigma))
  
      sigma = Symbol("sigma", real=True, positive=True)
      M = Moyal("M", mu, sigma)
      assert density(M)(z) == sqrt(2) * exp(
          -exp((mu - z) / sigma) / 2 - (-mu + z) / (2 * sigma)
      ) / (2 * sqrt(pi) * sigma)
      assert cdf(M)(z).simplify() == 1 - erf(sqrt(2) * exp((mu - z) / (2 * sigma)) / 2)
      assert characteristic_function(M)(z) == 2 ** (-I * sigma * z) * exp(
          I * mu * z
      ) * gamma(-I * sigma * z + Rational(1, 2)) / sqrt(pi)
      assert E(M) == mu + EulerGamma * sigma + sigma * log(2)
      assert moment_generating_function(M)(z) == 2 ** (-sigma * z) * exp(mu * z) * gamma(
          -sigma * z + Rational(1, 2)
      ) / sqrt(pi)
  
  
  def test_nakagami():
      mu = Symbol("mu", positive=True)
      omega = Symbol("omega", positive=True)
  
      X = Nakagami("x", mu, omega)
      assert density(X)(x) == (
          2
          * x ** (2 * mu - 1)
          * mu ** mu
          * omega ** (-mu)
          * exp(-(x ** 2) * mu / omega)
          / gamma(mu)
      )
      assert simplify(E(X)) == (
          sqrt(mu) * sqrt(omega) * gamma(mu + S.Half) / gamma(mu + 1)
      )
      assert simplify(variance(X)) == (
          omega - omega * gamma(mu + S.Half) ** 2 / (gamma(mu) * gamma(mu + 1))
      )
      assert cdf(X)(x) == Piecewise(
          (lowergamma(mu, mu * x ** 2 / omega) / gamma(mu), x > 0), (0, True)
      )
  
  
  def test_gaussian_inverse():
      # test for symbolic parameters
      a, b = symbols("a b")
      assert GaussianInverse("x", a, b)
  
      # Inverse Gaussian distribution is also known as Wald distribution
      # `GaussianInverse` can also be referred by the name `Wald`
      a, b, z = symbols("a b z")
      X = Wald("x", a, b)
      assert density(X)(z) == sqrt(2) * sqrt(b / z ** 3) * exp(
          -b * (-a + z) ** 2 / (2 * a ** 2 * z)
      ) / (2 * sqrt(pi))
  
      a, b = symbols("a b", positive=True)
      z = Symbol("z", positive=True)
  
      X = GaussianInverse("x", a, b)
      assert density(X)(z) == sqrt(2) * sqrt(b) * sqrt(z ** (-3)) * exp(
          -b * (-a + z) ** 2 / (2 * a ** 2 * z)
      ) / (2 * sqrt(pi))
      assert E(X) == a
      assert variance(X).expand() == a ** 3 / b
      assert (
          cdf(X)(z)
          == (S.Half - erf(sqrt(2) * sqrt(b) * (1 + z / a) / (2 * sqrt(z))) / 2)
          * exp(2 * b / a)
          + erf(sqrt(2) * sqrt(b) * (-1 + z / a) / (2 * sqrt(z))) / 2
          + S.Half
      )
  
      a = symbols("a", nonpositive=True)
      raises(ValueError, lambda: GaussianInverse("x", a, b))
  
      a = symbols("a", positive=True)
      b = symbols("b", nonpositive=True)
      raises(ValueError, lambda: GaussianInverse("x", a, b))
  
  
  def test_sampling_gaussian_inverse():
      scipy = import_module("scipy")
      if not scipy:
          skip("Scipy not installed. Abort tests for sampling of Gaussian inverse.")
      X = GaussianInverse("x", 1, 1)
      assert sample(X) in X.pspace.domain.set
  
  
  def test_pareto():
      xm, beta = symbols("xm beta", positive=True)
      alpha = beta + 5
      X = Pareto("x", xm, alpha)
  
      dens = density(X)
  
      # Tests cdf function
      assert cdf(X)(x) == Piecewise(
          (-(x ** (-beta - 5)) * xm ** (beta + 5) + 1, x >= xm), (0, True)
      )
  
      # Tests characteristic_function
      assert characteristic_function(X)(x) == (
          (-I * x * xm) ** (beta + 5) * (beta + 5) * uppergamma(-beta - 5, -I * x * xm)
      )
  
      assert dens(x) == x ** (-(alpha + 1)) * xm ** (alpha) * (alpha)
  
      assert simplify(E(X)) == alpha * xm / (alpha - 1)
  
      # computation of taylor series for MGF still too slow
      # assert simplify(variance(X)) == xm**2*alpha / ((alpha-1)**2*(alpha-2))
  
  
  def test_pareto_numeric():
      xm, beta = 3, 2
      alpha = beta + 5
      X = Pareto("x", xm, alpha)
  
      assert E(X) == alpha * xm / S(alpha - 1)
      assert variance(X) == xm ** 2 * alpha / S(((alpha - 1) ** 2 * (alpha - 2)))
      # Skewness tests too slow. Try shortcutting function?
  
  
  def test_PowerFunction():
      alpha = Symbol("alpha", nonpositive=True)
      a, b = symbols("a, b", real=True)
      raises(ValueError, lambda: PowerFunction("x", alpha, a, b))
  
      a, b = symbols("a, b", real=False)
      raises(ValueError, lambda: PowerFunction("x", alpha, a, b))
  
      alpha = Symbol("alpha", positive=True)
      a, b = symbols("a, b", real=True)
      raises(ValueError, lambda: PowerFunction("x", alpha, 5, 2))
  
      X = PowerFunction("X", 2, a, b)
      assert density(X)(z) == (-2 * a + 2 * z) / (-a + b) ** 2
      assert cdf(X)(z) == Piecewise(
          (
              a ** 2 / (a ** 2 - 2 * a * b + b ** 2)
              - 2 * a * z / (a ** 2 - 2 * a * b + b ** 2)
              + z ** 2 / (a ** 2 - 2 * a * b + b ** 2),
              a <= z,
          ),
          (0, True),
      )
  
      X = PowerFunction("X", 2, 0, 1)
      assert density(X)(z) == 2 * z
      assert cdf(X)(z) == Piecewise((z ** 2, z >= 0), (0, True))
      assert E(X) == Rational(2, 3)
      assert P(X < 0) == 0
      assert P(X < 1) == 1
  
  
  def test_raised_cosine():
      mu = Symbol("mu", real=True)
      s = Symbol("s", positive=True)
  
      X = RaisedCosine("x", mu, s)
  
      # Tests characteristics_function
      assert characteristic_function(X)(x) == Piecewise(
          (exp(-I * pi * mu / s) / 2, Eq(x, -pi / s)),
          (exp(I * pi * mu / s) / 2, Eq(x, pi / s)),
          (
              pi ** 2
              * exp(I * mu * x)
              * sin(s * x)
              / (s * x * (-(s ** 2) * x ** 2 + pi ** 2)),
              True,
          ),
      )
  
      assert density(X)(x) == (
          Piecewise(
              ((cos(pi * (x - mu) / s) + 1) / (2 * s), And(x <= mu + s, mu - s <= x)),
              (0, True),
          )
      )
  
  
  def test_rayleigh():
      sigma = Symbol("sigma", positive=True)
  
      X = Rayleigh("x", sigma)
  
      # Tests characteristic_function
      assert characteristic_function(X)(x) == (
          -sqrt(2)
          * sqrt(pi)
          * sigma
          * x
          * (erfi(sqrt(2) * sigma * x / 2) - I)
          * exp(-(sigma ** 2) * x ** 2 / 2)
          / 2
          + 1
      )
  
      assert density(X)(x) == x * exp(-(x ** 2) / (2 * sigma ** 2)) / sigma ** 2
      assert E(X) == sqrt(2) * sqrt(pi) * sigma / 2
      assert variance(X) == -pi * sigma ** 2 / 2 + 2 * sigma ** 2
      assert cdf(X)(x) == 1 - exp(-(x ** 2) / (2 * sigma ** 2))
      assert diff(cdf(X)(x), x) == density(X)(x)
  
  
  def test_reciprocal():
      a = Symbol("a", real=True)
      b = Symbol("b", real=True)
  
      X = Reciprocal("x", a, b)
      assert density(X)(x) == 1 / (x * (-log(a) + log(b)))
      assert cdf(X)(x) == Piecewise(
          (log(a) / (log(a) - log(b)) - log(x) / (log(a) - log(b)), a <= x), (0, True)
      )
      X = Reciprocal("x", 5, 30)
  
      assert E(X) == 25 / (log(30) - log(5))
      assert P(X < 4) == S.Zero
      assert P(X < 20) == log(20) / (log(30) - log(5)) - log(5) / (log(30) - log(5))
      assert cdf(X)(10) == log(10) / (log(30) - log(5)) - log(5) / (log(30) - log(5))
  
      a = symbols("a", nonpositive=True)
      raises(ValueError, lambda: Reciprocal("x", a, b))
  
      a = symbols("a", positive=True)
      b = symbols("b", positive=True)
      raises(ValueError, lambda: Reciprocal("x", a + b, a))
  
  
  def test_shiftedgompertz():
      b = Symbol("b", positive=True)
      eta = Symbol("eta", positive=True)
      X = ShiftedGompertz("x", b, eta)
      assert density(X)(x) == b * (eta * (1 - exp(-b * x)) + 1) * exp(-b * x) * exp(
          -eta * exp(-b * x)
      )
  
  
  def test_studentt():
      nu = Symbol("nu", positive=True)
  
      X = StudentT("x", nu)
      assert density(X)(x) == (1 + x ** 2 / nu) ** (-nu / 2 - S.Half) / (
          sqrt(nu) * beta(S.Half, nu / 2)
      )
      assert cdf(X)(x) == S.Half + x * gamma(nu / 2 + S.Half) * hyper(
          (S.Half, nu / 2 + S.Half), (Rational(3, 2),), -(x ** 2) / nu
      ) / (sqrt(pi) * sqrt(nu) * gamma(nu / 2))
  
  
  def test_trapezoidal():
      a = Symbol("a", real=True)
      b = Symbol("b", real=True)
      c = Symbol("c", real=True)
      d = Symbol("d", real=True)
  
      X = Trapezoidal("x", a, b, c, d)
      assert density(X)(x) == Piecewise(
          ((-2 * a + 2 * x) / ((-a + b) * (-a - b + c + d)), (a <= x) & (x < b)),
          (2 / (-a - b + c + d), (b <= x) & (x < c)),
          ((2 * d - 2 * x) / ((-c + d) * (-a - b + c + d)), (c <= x) & (x <= d)),
          (0, True),
      )
  
      X = Trapezoidal("x", 0, 1, 2, 3)
      assert E(X) == Rational(3, 2)
      assert variance(X) == Rational(5, 12)
      assert P(X < 2) == Rational(3, 4)
  
  
  def test_triangular():
      a = Symbol("a")
      b = Symbol("b")
      c = Symbol("c")
  
      X = Triangular("x", a, b, c)
      assert str(density(X)(x)) == (
          "Piecewise(((-2*a + 2*x)/((-a + b)*(-a + c)), (a <= x) & (c > x)), "
          "(2/(-a + b), Eq(c, x)), ((2*b - 2*x)/((-a + b)*(b - c)), (b >= x) & (c < x)), (0, True))"
      )
  
      # Tests moment_generating_function
      assert (
          moment_generating_function(X)(x).expand()
          == (
              (
                  -2 * (-a + b) * exp(c * x)
                  + 2 * (-a + c) * exp(b * x)
                  + 2 * (b - c) * exp(a * x)
              )
              / (x ** 2 * (-a + b) * (-a + c) * (b - c))
          ).expand()
      )
  
  
  def test_quadratic_u():
      a = Symbol("a", real=True)
      b = Symbol("b", real=True)
  
      X = QuadraticU("x", a, b)
      Y = QuadraticU("x", 1, 2)
  
      # Tests _moment_generating_function
      assert moment_generating_function(Y)(1) == -15 * exp(2) + 27 * exp(1)
      assert moment_generating_function(Y)(2) == -9 * exp(4) / 2 + 21 * exp(2) / 2
  
      assert density(X)(x) == (
          Piecewise(
              (12 * (x - a / 2 - b / 2) ** 2 / (-a + b) ** 3, And(x <= b, a <= x)),
              (0, True),
          )
      )
  
  
  def test_uniform():
      l = Symbol("l", real=True)
      w = Symbol("w", positive=True)
      X = Uniform("x", l, l + w)
  
      assert E(X) == l + w / 2
      assert variance(X).expand() == w ** 2 / 12
  
      # With numbers all is well
      X = Uniform("x", 3, 5)
      assert P(X < 3) == 0 and P(X > 5) == 0
      assert P(X < 4) == P(X > 4) == S.Half
  
      z = Symbol("z")
      p = density(X)(z)
      assert p.subs(z, 3.7) == S.Half
      assert p.subs(z, -1) == 0
      assert p.subs(z, 6) == 0
  
      c = cdf(X)
      assert c(2) == 0 and c(3) == 0
      assert c(Rational(7, 2)) == Rational(1, 4)
      assert c(5) == 1 and c(6) == 1
  
  
  @XFAIL
  def test_uniform_P():
      """ This stopped working because SingleContinuousPSpace.compute_density no
      longer calls integrate on a DiracDelta but rather just solves directly.
      integrate used to call UniformDistribution.expectation which special-cased
      subsed out the Min and Max terms that Uniform produces
  
      I decided to regress on this class for general cleanliness (and I suspect
      speed) of the algorithm.
      """
      l = Symbol("l", real=True)
      w = Symbol("w", positive=True)
      X = Uniform("x", l, l + w)
      assert P(X < l) == 0 and P(X > l + w) == 0
  
  
  def test_uniformsum():
      n = Symbol("n", integer=True)
      _k = Dummy("k")
      x = Symbol("x")
  
      X = UniformSum("x", n)
      res = Sum(
          (-1) ** _k * (-_k + x) ** (n - 1) * binomial(n, _k), (_k, 0, floor(x))
      ) / factorial(n - 1)
      assert density(X)(x).dummy_eq(res)
  
      # Tests set functions
      assert X.pspace.domain.set == Interval(0, n)
  
      # Tests the characteristic_function
      assert characteristic_function(X)(x) == (-I * (exp(I * x) - 1) / x) ** n
  
      # Tests the moment_generating_function
      assert moment_generating_function(X)(x) == ((exp(x) - 1) / x) ** n
  
  
  def test_von_mises():
      mu = Symbol("mu")
      k = Symbol("k", positive=True)
  
      X = VonMises("x", mu, k)
      assert density(X)(x) == exp(k * cos(x - mu)) / (2 * pi * besseli(0, k))
  
  
  def test_weibull():
      a, b = symbols("a b", positive=True)
      # FIXME: simplify(E(X)) seems to hang without extended_positive=True
      # On a Linux machine this had a rapid memory leak...
      # a, b = symbols('a b', positive=True)
      X = Weibull("x", a, b)
  
      assert E(X).expand() == a * gamma(1 + 1 / b)
      assert variance(X).expand() == (a ** 2 * gamma(1 + 2 / b) - E(X) ** 2).expand()
      assert simplify(skewness(X)) == (
          2 * gamma(1 + 1 / b) ** 3
          - 3 * gamma(1 + 1 / b) * gamma(1 + 2 / b)
          + gamma(1 + 3 / b)
      ) / (-gamma(1 + 1 / b) ** 2 + gamma(1 + 2 / b)) ** Rational(3, 2)
      assert (
          simplify(kurtosis(X))
          == (
              -3 * gamma(1 + 1 / b) ** 4
              + 6 * gamma(1 + 1 / b) ** 2 * gamma(1 + 2 / b)
              - 4 * gamma(1 + 1 / b) * gamma(1 + 3 / b)
              + gamma(1 + 4 / b)
          )
          / (gamma(1 + 1 / b) ** 2 - gamma(1 + 2 / b)) ** 2
      )
  
  
  def test_weibull_numeric():
      # Test for integers and rationals
      a = 1
      bvals = [S.Half, 1, Rational(3, 2), 5]
      for b in bvals:
          X = Weibull("x", a, b)
          assert simplify(E(X)) == expand_func(a * gamma(1 + 1 / S(b)))
          assert simplify(variance(X)) == simplify(
              a ** 2 * gamma(1 + 2 / S(b)) - E(X) ** 2
          )
          # Not testing Skew... it's slow with int/frac values > 3/2
  
  
  def test_wignersemicircle():
      R = Symbol("R", positive=True)
  
      X = WignerSemicircle("x", R)
      assert density(X)(x) == 2 * sqrt(-(x ** 2) + R ** 2) / (pi * R ** 2)
      assert E(X) == 0
  
      # Tests ChiNoncentralDistribution
      assert characteristic_function(X)(x) == Piecewise(
          (2 * besselj(1, R * x) / (R * x), Ne(x, 0)), (1, True)
      )
  
  
  def test_prefab_sampling():
      N = Normal("X", 0, 1)
      L = LogNormal("L", 0, 1)
      E = Exponential("Ex", 1)
      P = Pareto("P", 1, 3)
      W = Weibull("W", 1, 1)
      U = Uniform("U", 0, 1)
      B = Beta("B", 2, 5)
      G = Gamma("G", 1, 3)
  
      variables = [N, L, E, P, W, U, B, G]
      niter = 10
      for var in variables:
          for i in range(niter):
              assert sample(var) in var.pspace.domain.set
  
  
  def test_input_value_assertions():
      a, b = symbols("a b")
      p, q = symbols("p q", positive=True)
      m, n = symbols("m n", positive=False, real=True)
  
      raises(ValueError, lambda: Normal("x", 3, 0))
      raises(ValueError, lambda: Normal("x", m, n))
      Normal("X", a, p)  # No error raised
      raises(ValueError, lambda: Exponential("x", m))
      Exponential("Ex", p)  # No error raised
      for fn in [Pareto, Weibull, Beta, Gamma]:
          raises(ValueError, lambda: fn("x", m, p))
          raises(ValueError, lambda: fn("x", p, n))
          fn("x", p, q)  # No error raised
  
  
  def test_unevaluated():
      X = Normal("x", 0, 1)
      assert str(E(X, evaluate=False)) == (
          "Integral(sqrt(2)*x*exp(-x**2/2)/" "(2*sqrt(pi)), (x, -oo, oo))"
      )
  
      assert str(E(X + 1, evaluate=False)) == (
          "Integral(sqrt(2)*x*exp(-x**2/2)/" "(2*sqrt(pi)), (x, -oo, oo)) + 1"
      )
  
      assert str(P(X > 0, evaluate=False)) == (
          "Integral(sqrt(2)*exp(-_z**2/2)/" "(2*sqrt(pi)), (_z, 0, oo))"
      )
  
      assert P(X > 0, X ** 2 < 1, evaluate=False) == S.Half
  
  
  def test_probability_unevaluated():
      T = Normal("T", 30, 3)
      assert type(P(T > 33, evaluate=False)) == Integral
  
  
  def test_density_unevaluated():
      X = Normal("X", 0, 1)
      Y = Normal("Y", 0, 2)
      assert isinstance(density(X + Y, evaluate=False)(z), Integral)
  
  
  def test_NormalDistribution():
      nd = NormalDistribution(0, 1)
      x = Symbol("x")
      assert nd.cdf(x) == erf(sqrt(2) * x / 2) / 2 + S.Half
      assert isinstance(nd.sample(), float) or nd.sample().is_Number
      assert nd.expectation(1, x) == 1
      assert nd.expectation(x, x) == 0
      assert nd.expectation(x ** 2, x) == 1
      # Test issue 10076
      a = SingleContinuousPSpace(x, NormalDistribution(2, 4))
      _z = Dummy("_z")
  
      expected1 = Integral(
          sqrt(2) * exp(-((_z - 2) ** 2) / 32) / (8 * sqrt(pi)), (_z, -oo, 1)
      )
      assert a.probability(x < 1, evaluate=False).dummy_eq(expected1) is True
  
      expected2 = Integral(
          sqrt(2) * exp(-((_z - 2) ** 2) / 32) / (8 * sqrt(pi)), (_z, 1, oo)
      )
      assert a.probability(x > 1, evaluate=False).dummy_eq(expected2) is True
  
      b = SingleContinuousPSpace(x, NormalDistribution(1, 9))
  
      expected3 = Integral(
          sqrt(2) * exp(-((_z - 1) ** 2) / 162) / (18 * sqrt(pi)), (_z, 6, oo)
      )
      assert b.probability(x > 6, evaluate=False).dummy_eq(expected3) is True
  
      expected4 = Integral(
          sqrt(2) * exp(-((_z - 1) ** 2) / 162) / (18 * sqrt(pi)), (_z, -oo, 6)
      )
      assert b.probability(x < 6, evaluate=False).dummy_eq(expected4) is True
  
  
  def test_random_parameters():
      mu = Normal("mu", 2, 3)
      meas = Normal("T", mu, 1)
      assert density(meas, evaluate=False)(z)
      assert isinstance(pspace(meas), JointPSpace)
      # assert density(meas, evaluate=False)(z) == Integral(mu.pspace.pdf *
      #        meas.pspace.pdf, (mu.symbol, -oo, oo)).subs(meas.symbol, z)
  
  
  def test_random_parameters_given():
      mu = Normal("mu", 2, 3)
      meas = Normal("T", mu, 1)
      assert given(meas, Eq(mu, 5)) == Normal("T", 5, 1)
  
  
  def test_conjugate_priors():
      mu = Normal("mu", 2, 3)
      x = Normal("x", mu, 1)
      assert isinstance(simplify(density(mu, Eq(x, y), evaluate=False)(z)), Mul)
  
  
  def test_difficult_univariate():
      """ Since using solve in place of deltaintegrate we're able to perform
      substantially more complex density computations on single continuous random
      variables """
      x = Normal("x", 0, 1)
      assert density(x ** 3)
      assert density(exp(x ** 2))
      assert density(log(x))
  
  
  def test_issue_10003():
      X = Exponential("x", 3)
      G = Gamma("g", 1, 2)
      assert P(X < -1) is S.Zero
      assert P(G < -1) is S.Zero
  
  
  @slow
  def test_precomputed_cdf():
      x = symbols("x", real=True)
      mu = symbols("mu", real=True)
      sigma, xm, alpha = symbols("sigma xm alpha", positive=True)
      n = symbols("n", integer=True, positive=True)
      distribs = [
          Normal("X", mu, sigma),
          Pareto("P", xm, alpha),
          ChiSquared("C", n),
          Exponential("E", sigma),
          # LogNormal("L", mu, sigma),
      ]
      for X in distribs:
          compdiff = cdf(X)(x) - simplify(X.pspace.density.compute_cdf()(x))
          compdiff = simplify(compdiff.rewrite(erfc))
          assert compdiff == 0
  
  
  @slow
  def test_precomputed_characteristic_functions():
      import mpmath
  
      def test_cf(dist, support_lower_limit, support_upper_limit):
          pdf = density(dist)
          t = Symbol("t")
  
          # first function is the hardcoded CF of the distribution
          cf1 = lambdify([t], characteristic_function(dist)(t), "mpmath")
  
          # second function is the Fourier transform of the density function
          f = lambdify([x, t], pdf(x) * exp(I * x * t), "mpmath")
          cf2 = lambda t: mpmath.quad(
              lambda x: f(x, t), [support_lower_limit, support_upper_limit], maxdegree=10
          )
  
          # compare the two functions at various points
          for test_point in [2, 5, 8, 11]:
              n1 = cf1(test_point)
              n2 = cf2(test_point)
  
              assert abs(re(n1) - re(n2)) < 1e-12
              assert abs(im(n1) - im(n2)) < 1e-12
  
      test_cf(Beta("b", 1, 2), 0, 1)
      test_cf(Chi("c", 3), 0, mpmath.inf)
      test_cf(ChiSquared("c", 2), 0, mpmath.inf)
      test_cf(Exponential("e", 6), 0, mpmath.inf)
      test_cf(Logistic("l", 1, 2), -mpmath.inf, mpmath.inf)
      test_cf(Normal("n", -1, 5), -mpmath.inf, mpmath.inf)
      test_cf(RaisedCosine("r", 3, 1), 2, 4)
      test_cf(Rayleigh("r", 0.5), 0, mpmath.inf)
      test_cf(Uniform("u", -1, 1), -1, 1)
      test_cf(WignerSemicircle("w", 3), -3, 3)
  
  
  def test_long_precomputed_cdf():
      x = symbols("x", real=True)
      distribs = [
          Arcsin("A", -5, 9),
          Dagum("D", 4, 10, 3),
          Erlang("E", 14, 5),
          Frechet("F", 2, 6, -3),
          Gamma("G", 2, 7),
          GammaInverse("GI", 3, 5),
          Kumaraswamy("K", 6, 8),
          Laplace("LA", -5, 4),
          Logistic("L", -6, 7),
          Nakagami("N", 2, 7),
          StudentT("S", 4),
      ]
      for distr in distribs:
          for _ in range(5):
              assert tn(diff(cdf(distr)(x), x), density(distr)(x), x, a=0, b=0, c=1, d=0)
  
      US = UniformSum("US", 5)
      pdf01 = density(US)(x).subs(floor(x), 0).doit()  # pdf on (0, 1)
      cdf01 = cdf(US, evaluate=False)(x).subs(floor(x), 0).doit()  # cdf on (0, 1)
      assert tn(diff(cdf01, x), pdf01, x, a=0, b=0, c=1, d=0)
  
  
  def test_issue_13324():
      X = Uniform("X", 0, 1)
      assert E(X, X > S.Half) == Rational(3, 4)
      assert E(X, X > 0) == S.Half
  
  
  def test_FiniteSet_prob():
      E = Exponential("E", 3)
      N = Normal("N", 5, 7)
      assert P(Eq(E, 1)) is S.Zero
      assert P(Eq(N, 2)) is S.Zero
      assert P(Eq(N, x)) is S.Zero
  
  
  def test_prob_neq():
      E = Exponential("E", 4)
      X = ChiSquared("X", 4)
      assert P(Ne(E, 2)) == 1
      assert P(Ne(X, 4)) == 1
      assert P(Ne(X, 4)) == 1
      assert P(Ne(X, 5)) == 1
      assert P(Ne(E, x)) == 1
  
  
  def test_union():
      N = Normal("N", 3, 2)
      assert simplify(P(N ** 2 - N > 2)) == -erf(sqrt(2)) / 2 - erfc(
          sqrt(2) / 4
      ) / 2 + Rational(3, 2)
      assert simplify(P(N ** 2 - 4 > 0)) == -erf(5 * sqrt(2) / 4) / 2 - erfc(
          sqrt(2) / 4
      ) / 2 + Rational(3, 2)
  
  
  def test_Or():
      N = Normal("N", 0, 1)
      assert simplify(P(Or(N > 2, N < 1))) == -erf(sqrt(2)) / 2 - erfc(
          sqrt(2) / 2
      ) / 2 + Rational(3, 2)
      assert P(Or(N < 0, N < 1)) == P(N < 1)
      assert P(Or(N > 0, N < 0)) == 1
  
  
  def test_conditional_eq():
      E = Exponential("E", 1)
      assert P(Eq(E, 1), Eq(E, 1)) == 1
      assert P(Eq(E, 1), Eq(E, 2)) == 0
      assert P(E > 1, Eq(E, 2)) == 1
      assert P(E < 1, Eq(E, 2)) == 0
  