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)))

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_beta():
    a, b, x = symbols("a b x", positive=True)
    e = x**(a - 1)*(-x + 1)**(b - 1)/beta(a, b)
    Q = QQ[a, b].get_field()
    h1 = expr_to_holonomic(e, x, domain=Q)

    _, Dx = DifferentialOperators(Q.old_poly_ring(x), 'Dx')
    h2 = HolonomicFunction((a + x*(-a - b + 2) - 1) + (x**2 - x)*Dx, x)

    assert h1 == h2

def test_hyper():
    for x in sorted(exparg):
        test("erf", x, N(sp.erf(x)))
    for x in sorted(exparg):
        test("erfc", x, N(sp.erfc(x)))

    gamarg = FiniteSet(*(x+S(1)/12 for x in exparg))
    betarg = ProductSet(gamarg, gamarg)
    for x in sorted(gamarg):
        test("lgamma", x, N(sp.log(abs(sp.gamma(x)))))
    for x in sorted(gamarg):
        test("gamma", x, N(sp.gamma(x)))
    for x, y in sorted(betarg, key=lambda (x, y): (y, x)):
        test("beta", x, y, N(sp.beta(x, y)))

    pgamarg = FiniteSet(S(1)/12, S(1)/3, S(3)/2, 5)
    pgamargp = ProductSet(gamarg & Interval(0, oo, True), pgamarg)
    for a, x in sorted(pgamargp):
        test("pgamma", a, x, N(sp.lowergamma(a, x)))
    for a, x in sorted(pgamargp):
        test("pgammac", a, x, N(sp.uppergamma(a, x)))
    for a, x in sorted(pgamargp):
        test("pgammar", a, x, N(sp.lowergamma(a, x)/sp.gamma(a)))
    for a, x in sorted(pgamargp):
        test("pgammarc", a, x, N(sp.uppergamma(a, x)/sp.gamma(a)))
    for a, x in sorted(pgamargp):
        test("ipgammarc", a, N(sp.uppergamma(a, x)/sp.gamma(a)), x)

    pbetargp = [(a, b, x) for a, b, x in ProductSet(betarg, pgamarg)
                if a > 0 and b > 0 and x < 1]
    pbetargp.sort(key=lambda (a, b, x): (b, a, x))
    for a, b, x in pbetargp:
        test("pbeta", a, b, x, mp.betainc(mpf(a), mpf(b), x2=mpf(x)))
    for a, b, x in pbetargp:
        test("pbetar", a, b, x, mp.betainc(mpf(a), mpf(b), x2=mpf(x),
                                           regularized=True))
    for a, b, x in pbetargp:
        test("ipbetar", a, b, mp.betainc(mpf(a), mpf(b), x2=mpf(x),
                                         regularized=True), x)

    for x in sorted(posarg):
        test("j0", x, N(sp.besselj(0, x)))
    for x in sorted(posarg):
        test("j1", x, N(sp.besselj(1, x)))
    for x in sorted(posarg-FiniteSet(0)):
        test("y0", x, N(sp.bessely(0, x)))
    for x in sorted(posarg-FiniteSet(0)):
        test("y1", x, N(sp.bessely(1, x)))

def test_beta():
    a, b = symbols('alpha beta', positive=True)

    B = Beta(a, b)

    assert pspace(B).domain.set == Interval(0, 1)

    x, dens = Density(B)
    assert dens == x**(a-1)*(1-x)**(b-1) / beta(a,b)

    # This is too slow
    # assert E(B) == a / (a + b)
    # assert Var(B) == (a*b) / ((a+b)**2 * (a+b+1))

    # Full symbolic solution is too much, test with numeric version
    a, b = 1, 2
    B = Beta(a, b)
    assert E(B) == a / S(a + b)
    assert Var(B) == (a*b) / S((a+b)**2 * (a+b+1))

def test_studentt():
    nu = Symbol("nu", positive=True)

    X = StudentT('x', nu)
    assert density(X)(x) == (1 + x**2/nu)**(-nu/2 - 1/2)/(sqrt(nu)*beta(1/2, nu/2))

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)

def test_beta_function():
    x, y = Symbol('x'), Symbol('y')
    assert beta(x,y) == gamma(x)*gamma(y)/gamma(x+y)
    assert beta(x,y) == beta(y,x) # Symmetric

 def _cdf(self, x):
     return Piecewise((1 - floor(x) * beta(floor(x), self.rho + 1), x >= 1), (0, True))

 def pdf(self, k):
     rho = self.rho
     return rho * beta(k, rho + 1)

def test_beta():
    x, y = Symbol('x'), Symbol('y')

    assert isinstance(beta(x, y), beta)

    assert expand_func(beta(x, y)) == gamma(x)*gamma(y)/gamma(x + y)
    assert expand_func(beta(x, y) - beta(y, x)) == 0  # Symmetric
    assert expand_func(beta(x, y)) == expand_func(beta(x, y + 1) + beta(x + 1, y)).simplify()

    assert diff(beta(x, y), x) == beta(x, y)*(digamma(x) - digamma(x + y))
    assert diff(beta(x, y), y) == beta(x, y)*(digamma(y) - digamma(x + y))