Python scipy.special.hyp2f1() Examples

def test_hyp2f1_some_points_2():
    # Taken from mpmath unit tests -- this point failed for mpmath 0.13 but
    # was fixed in their SVN since then
    pts = [
        (112, (51,10), (-9,10), -0.99999),
        (10,-900,10.5,0.99),
        (10,-900,-10.5,0.99),
    ]

    def fev(x):
        if isinstance(x, tuple):
            return float(x[0]) / x[1]
        else:
            return x

    dataset = [tuple(map(fev, p)) + (float(mpmath.hyp2f1(*p)),) for p in pts]
    dataset = np.array(dataset, dtype=np.float_)

    FuncData(sc.hyp2f1, dataset, (0,1,2,3), 4, rtol=1e-10).check() 

def test_hyp2f1_real_some():
    dataset = []
    for a in [-10, -5, -1.8, 1.8, 5, 10]:
        for b in [-2.5, -1, 1, 7.4]:
            for c in [-9, -1.8, 5, 20.4]:
                for z in [-10, -1.01, -0.99, 0, 0.6, 0.95, 1.5, 10]:
                    try:
                        v = float(mpmath.hyp2f1(a, b, c, z))
                    except:
                        continue
                    dataset.append((a, b, c, z, v))
    dataset = np.array(dataset, dtype=np.float_)

    olderr = np.seterr(invalid='ignore')
    try:
        FuncData(sc.hyp2f1, dataset, (0,1,2,3), 4, rtol=1e-9,
                 ignore_inf_sign=True).check()
    finally:
        np.seterr(**olderr) 

def derivatives(self, x, y, b, t, q):
        """
        returns the deflection
        """
        # elliptical radius, eq. (5)
        Z = np.empty(np.shape(x), dtype=complex)
        Z.real = q*x
        Z.imag = y
        R = np.abs(Z)

        # angular dependency with extra factor of R, eq. (23)
        R_omega = Z*hyp2f1(1, t/2, 2-t/2, -(1-q)/(1+q)*(Z/Z.conj()))

        # deflection, eq. (22)
        alpha = 2/(1+q)*(b/R)**t*R_omega

        # return real and imaginary part
        return alpha.real, alpha.imag 

def expected_number_of_purchases_up_to_time(self, t):
        """
        Return expected number of repeat purchases up to time t.

        Calculate the expected number of repeat purchases up to time t for a
        randomly choose individual from the population.

        Parameters
        ----------
        t: array_like
            times to calculate the expectation for

        Returns
        -------
        array_like

        """
        r, alpha, a, b = self._unload_params("r", "alpha", "a", "b")
        hyp = hyp2f1(r, b + 1, a + b, t / (alpha + t))
        return b / (a - 1) * (1 - hyp * (alpha / (alpha + t)) ** r) 

def test_hyp2f1_real_some():
    dataset = []
    for a in [-10, -5, -1.8, 1.8, 5, 10]:
        for b in [-2.5, -1, 1, 7.4]:
            for c in [-9, -1.8, 5, 20.4]:
                for z in [-10, -1.01, -0.99, 0, 0.6, 0.95, 1.5, 10]:
                    try:
                        v = float(mpmath.hyp2f1(a, b, c, z))
                    except:
                        continue
                    dataset.append((a, b, c, z, v))
    dataset = np.array(dataset, dtype=np.float_)

    olderr = np.seterr(invalid='ignore')
    try:
        FuncData(sc.hyp2f1, dataset, (0,1,2,3), 4, rtol=1e-9,
                 ignore_inf_sign=True).check()
    finally:
        np.seterr(**olderr) 

def expected_number_of_purchases_up_to_time(
        self, 
        t
    ):
        """
        Calculate the expected number of repeat purchases up to time t.

        Calculate repeat purchases for a randomly chosen individual from the
        population.

        Equivalent to equation (9) of [2]_.

        Parameters
        ----------
        t: array_like
            times to calculate the expection for

        Returns
        -------
        array_like

        References
        ----------
        .. [2] Fader, Peter S., Bruce G.S. Hardie, and Ka Lok Lee (2005a),
        "Counting Your Customers the Easy Way: An Alternative to the
        Pareto/NBD Model," Marketing Science, 24 (2), 275-84.
        """

        r, alpha, a, b = self._unload_params("r", "alpha", "a", "b")
        hyp = hyp2f1(r, b, a + b - 1, t / (alpha + t))

        return (a + b - 1) / (a - 1) * (1 - hyp * (alpha / (alpha + t)) ** r) 

def test_hyp2f1_real_random():
    npoints = 500
    dataset = np.zeros((npoints, 5), np.float_)

    np.random.seed(1234)
    dataset[:, 0] = np.random.pareto(1.5, npoints)
    dataset[:, 1] = np.random.pareto(1.5, npoints)
    dataset[:, 2] = np.random.pareto(1.5, npoints)
    dataset[:, 3] = 2*np.random.rand(npoints) - 1

    dataset[:, 0] *= (-1)**np.random.randint(2, npoints)
    dataset[:, 1] *= (-1)**np.random.randint(2, npoints)
    dataset[:, 2] *= (-1)**np.random.randint(2, npoints)

    for ds in dataset:
        if mpmath.__version__ < '0.14':
            # mpmath < 0.14 fails for c too much smaller than a, b
            if abs(ds[:2]).max() > abs(ds[2]):
                ds[2] = abs(ds[:2]).max()
        ds[4] = float(mpmath.hyp2f1(*tuple(ds[:4])))

    FuncData(sc.hyp2f1, dataset, (0, 1, 2, 3), 4, rtol=1e-9).check()

# ------------------------------------------------------------------------------
# erf (complex)
# ------------------------------------------------------------------------------ 

def test_hyp2f1_complex(self):
        # Scipy's hyp2f1 seems to have performance and accuracy problems
        assert_mpmath_equal(lambda a, b, c, x: sc.hyp2f1(a.real, b.real, c.real, x),
                            exception_to_nan(lambda a, b, c, x: mpmath.hyp2f1(a, b, c, x, **HYPERKW)),
                            [Arg(-1e2, 1e2), Arg(-1e2, 1e2), Arg(-1e2, 1e2), ComplexArg()],
                            n=10) 

def tpl_exp_spec_dens(k, dim, len_scale, hurst, len_low=0.0):
    r"""
    Spectal density of the TPLExponential covariance model.

    Parameters
    ----------
    k : :class:`float`
        Radius of the phase: :math:`k=\left\Vert\mathbf{k}\right\Vert`
    dim : :class:`int`
        Dimension of the model.
    len_scale : :class:`float`
        Length scale of the model.
    hurst : :class:`float`
        Hurst coefficient of the power law.
    len_low : :class:`float`, optional
        The lower length scale truncation of the model.
        Default: 0.0

    Returns
    -------
    :class:`float`
        spectal density of the TPLExponential model
    """
    if np.isclose(len_low, 0.0):
        k = np.array(k, dtype=np.double)
        z = (k * len_scale) ** 2
        a = hurst + dim / 2.0
        b = hurst + 0.5
        c = hurst + dim / 2.0 + 1.0
        d = dim / 2.0 + 0.5
        fac = len_scale ** dim * hurst * sps.gamma(d) / (np.pi ** d * a)
        return fac / (1.0 + z) ** a * sps.hyp2f1(a, b, c, z / (1.0 + z))
    fac_up = (len_scale + len_low) ** (2 * hurst)
    spec_up = tpl_exp_spec_dens(k, dim, len_scale + len_low, hurst)
    fac_low = len_low ** (2 * hurst)
    spec_low = tpl_exp_spec_dens(k, dim, len_low, hurst)
    return (fac_up * spec_up - fac_low * spec_low) / (fac_up - fac_low) 

def hyp2f1(*args, **kwargs):
        from scipy.special import hyp2f1
        return hyp2f1(*args, **kwargs) 

def erf(*args, **kwargs):
            from scipy.special import erf
            return erf(*args, **kwargs)


#    from scipy.special import lambertw, ellipe, gammaincc, gamma # fluids
#    from scipy.special import i1, i0, k1, k0, iv # ht
#    from scipy.special import hyp2f1    
#    if erf is None:
#        from scipy.special import erf 

def hyp2f1(*args, **kwargs):
        import mpmath
        return mpmath.hyp2f1(*args, **kwargs) 

def _log_A_0(
        params, 
        freq, 
        recency, 
        age
    ):
        """
        log_A_0.
        
        Equation (19) and (20) from paper:
        http://brucehardie.com/notes/009/pareto_nbd_derivations_2005-11-05.pdf
        """

        r, alpha, s, beta = params

        if alpha < beta:
            min_of_alpha_beta, max_of_alpha_beta, t = (alpha, beta, r + freq)
        else:
            min_of_alpha_beta, max_of_alpha_beta, t = (beta, alpha, s + 1)
        abs_alpha_beta = max_of_alpha_beta - min_of_alpha_beta

        rsf = r + s + freq
        p_1 = hyp2f1(rsf, t, rsf + 1.0, abs_alpha_beta / (max_of_alpha_beta + recency))
        q_1 = max_of_alpha_beta + recency
        p_2 = hyp2f1(rsf, t, rsf + 1.0, abs_alpha_beta / (max_of_alpha_beta + age))
        q_2 = max_of_alpha_beta + age

        try:
            size = len(freq)
            sign = np.ones(size)
        except TypeError:
            sign = 1

        return logsumexp([log(p_1) + rsf * log(q_2), log(p_2) + rsf * log(q_1)], axis=0, b=[sign, -sign]) - rsf * log(
            q_1 * q_2
        ) 

def _T_hypergeometric(self, x):
        """ Compute T_hypergeometric(x) using Gauss Hypergeometric function 2F1

        T(x) = 2 \\sqrt(x) _{2}F_{1} \\left(\\frac{1}{6}, \\frac{1}{2}; \\frac{7}{6}; -x^3)

        Note:
        The scipy.special.hyp2f1 code already implements the hypergeometric
        transformation suggested by
            Baes, Camps, Van De Putte, 2017, MNRAS, 468, 927.
        for use in actual numerical evaulations.

        """
        from scipy.special import hyp2f1
        return 2 * np.sqrt(x) * hyp2f1(1./6, 1./2, 7./6, -x**3) 

def conditional_expected_number_of_purchases_up_to_time(self, t, frequency, recency, T):
        """
        Conditional expected number of repeat purchases up to time t.

        Calculate the expected number of repeat purchases up to time t for a
        randomly choose individual from the population, given they have
        purchase history (frequency, recency, T)
        See Wagner, U. and Hoppe D. (2008).

        Parameters
        ----------
        t: array_like
            times to calculate the expectation for.
        frequency: array_like
            historical frequency of customer.
        recency: array_like
            historical recency of customer.
        T: array_like
            age of the customer.

        Returns
        -------
        array_like

        """
        x = frequency
        r, alpha, a, b = self._unload_params("r", "alpha", "a", "b")

        hyp_term = hyp2f1(r + x, b + x + 1, a + b + x, t / (alpha + T + t))
        first_term = (a + b + x) / (a - 1)
        second_term = 1 - hyp_term * ((alpha + T) / (alpha + t + T)) ** (r + x)
        numerator = first_term * second_term

        denominator = 1 + (a / (b + x)) * ((alpha + T) / (alpha + recency)) ** (r + x)

        return numerator / denominator 

def resid_anscombe(self, endog, mu):
        """
        The Anscombe residuals for the negative binomial family

        Parameters
        ----------
        endog : array-like
            Endogenous response variable
        mu : array-like
            Fitted mean response variable

        Returns
        -------
        resid_anscombe : array
            The Anscombe residuals as defined below.

        Notes
        -----
        `resid_anscombe` = (hyp2f1(-alpha*endog)-hyp2f1(-alpha*mu)+\
                1.5*(endog**(2/3.)-mu**(2/3.)))/(mu+alpha*mu**2)**(1/6.)

        where hyp2f1 is the hypergeometric 2f1 function parameterized as
        hyp2f1(x) = hyp2f1(2/3.,1/3.,5/3.,x)
        """

        hyp2f1 = lambda x: special.hyp2f1(2 / 3., 1 / 3., 5 / 3., x)
        return ((hyp2f1(-self.alpha * endog) - hyp2f1(-self.alpha * mu) +
                 1.5 * (endog**(2 / 3.) - mu**(2 / 3.))) /
                (mu + self.alpha * mu**2)**(1 / 6.)) 

def _pdf(self, x, a, b, c, z):
        Cinv = sc.gamma(a)*sc.gamma(b)/sc.gamma(a+b)*sc.hyp2f1(c, a, a+b, -z)
        return 1.0/Cinv * x**(a-1.0) * (1.0-x)**(b-1.0) / (1.0+z*x)**c 

def _munp(self, n, a, b, c, z):
        fac = sc.beta(n+a, b) / sc.beta(a, b)
        num = sc.hyp2f1(c, a+n, a+b+n, -z)
        den = sc.hyp2f1(c, a, a+b, -z)
        return fac*num / den 

def _cdf(self, x, c):
        term1 = x / sc.beta(0.5, c / 2.0)
        res = 0.5 + term1 * sc.hyp2f1(0.5, 1 - c / 2.0, 1.5, x**2)
        # There's an issue with hyp2f1, it returns nans near x = +-1, c > 100.
        # Use the generic implementation in that case.  See gh-1285 for
        # background.
        if np.any(np.isnan(res)):
            return rv_continuous._cdf(self, x, c)
        return res 

def _pdf(self, x, a, b, c, z):
        # gausshyper.pdf(x, a, b, c, z) =
        #   C * x**(a-1) * (1-x)**(b-1) * (1+z*x)**(-c)
        Cinv = sc.gamma(a)*sc.gamma(b)/sc.gamma(a+b)*sc.hyp2f1(c, a, a+b, -z)
        return 1.0/Cinv * x**(a-1.0) * (1.0-x)**(b-1.0) / (1.0+z*x)**c 

def test_hyp2f1_real_random():
    dataset = []

    npoints = 500
    dataset = np.zeros((npoints, 5), np.float_)

    np.random.seed(1234)
    dataset[:,0] = np.random.pareto(1.5, npoints)
    dataset[:,1] = np.random.pareto(1.5, npoints)
    dataset[:,2] = np.random.pareto(1.5, npoints)
    dataset[:,3] = 2*np.random.rand(npoints) - 1

    dataset[:,0] *= (-1)**np.random.randint(2, npoints)
    dataset[:,1] *= (-1)**np.random.randint(2, npoints)
    dataset[:,2] *= (-1)**np.random.randint(2, npoints)

    for ds in dataset:
        if mpmath.__version__ < '0.14':
            # mpmath < 0.14 fails for c too much smaller than a, b
            if abs(ds[:2]).max() > abs(ds[2]):
                ds[2] = abs(ds[:2]).max()
        ds[4] = float(mpmath.hyp2f1(*tuple(ds[:4])))

    FuncData(sc.hyp2f1, dataset, (0,1,2,3), 4, rtol=1e-9).check()

#------------------------------------------------------------------------------
# erf (complex)
#------------------------------------------------------------------------------ 

def _munp(self, n, a, b, c, z):
        fac = sc.beta(n+a, b) / sc.beta(a, b)
        num = sc.hyp2f1(c, a+n, a+b+n, -z)
        den = sc.hyp2f1(c, a, a+b, -z)
        return fac*num / den 

def _cdf(self, x, c):
        term1 = x / sc.beta(0.5, c / 2.0)
        res = 0.5 + term1 * sc.hyp2f1(0.5, 1 - c / 2.0, 1.5, x**2)
        # There's an issue with hyp2f1, it returns nans near x = +-1, c > 100.
        # Use the generic implementation in that case.  See gh-1285 for
        # background.
        if np.any(np.isnan(res)):
            return rv_continuous._cdf(self, x, c)
        return res 

def test_hyp2f1(self):
        assert_equal(cephes.hyp2f1(1,1,1,0),1.0) 

def test_hyp2f1(self):
        # a collection of special cases taken from AMS 55
        values = [[0.5, 1, 1.5, 0.2**2, 0.5/0.2*log((1+0.2)/(1-0.2))],
                  [0.5, 1, 1.5, -0.2**2, 1./0.2*arctan(0.2)],
                  [1, 1, 2, 0.2, -1/0.2*log(1-0.2)],
                  [3, 3.5, 1.5, 0.2**2,
                      0.5/0.2/(-5)*((1+0.2)**(-5)-(1-0.2)**(-5))],
                  [-3, 3, 0.5, sin(0.2)**2, cos(2*3*0.2)],
                  [3, 4, 8, 1, special.gamma(8)*special.gamma(8-4-3)/special.gamma(8-3)/special.gamma(8-4)],
                  [3, 2, 3-2+1, -1, 1./2**3*sqrt(pi) *
                      special.gamma(1+3-2)/special.gamma(1+0.5*3-2)/special.gamma(0.5+0.5*3)],
                  [5, 2, 5-2+1, -1, 1./2**5*sqrt(pi) *
                      special.gamma(1+5-2)/special.gamma(1+0.5*5-2)/special.gamma(0.5+0.5*5)],
                  [4, 0.5+4, 1.5-2*4, -1./3, (8./9)**(-2*4)*special.gamma(4./3) *
                      special.gamma(1.5-2*4)/special.gamma(3./2)/special.gamma(4./3-2*4)],
                  # and some others
                  # ticket #424
                  [1.5, -0.5, 1.0, -10.0, 4.1300097765277476484],
                  # negative integer a or b, with c-a-b integer and x > 0.9
                  [-2,3,1,0.95,0.715],
                  [2,-3,1,0.95,-0.007],
                  [-6,3,1,0.95,0.0000810625],
                  [2,-5,1,0.95,-0.000029375],
                  # huge negative integers
                  (10, -900, 10.5, 0.99, 1.91853705796607664803709475658e-24),
                  (10, -900, -10.5, 0.99, 3.54279200040355710199058559155e-18),
                  ]
        for i, (a, b, c, x, v) in enumerate(values):
            cv = special.hyp2f1(a, b, c, x)
            assert_almost_equal(cv, v, 8, err_msg='test #%d' % i) 

def test_expi_complex():
    dataset = []
    for r in np.logspace(-99, 2, 10):
        for p in np.linspace(0, 2*np.pi, 30):
            z = r*np.exp(1j*p)
            dataset.append((z, complex(mpmath.ei(z))))
    dataset = np.array(dataset, dtype=np.complex_)

    FuncData(sc.expi, dataset, 0, 1).check()


#------------------------------------------------------------------------------
# hyp2f1
#------------------------------------------------------------------------------ 

def test_hyp2f1_strange_points():
    pts = [
        (2,-1,-1,0.7),
        (2,-2,-2,0.7),
    ]
    kw = dict(eliminate=True)
    dataset = [p + (float(mpmath.hyp2f1(*p, **kw)),) for p in pts]
    dataset = np.array(dataset, dtype=np.float_)

    FuncData(sc.hyp2f1, dataset, (0,1,2,3), 4, rtol=1e-10).check() 

def test_hyp2f1_real_some_points():
    pts = [
        (1,2,3,0),
        (1./3, 2./3, 5./6, 27./32),
        (1./4, 1./2, 3./4, 80./81),
        (2,-2,-3,3),
        (2,-3,-2,3),
        (2,-1.5,-1.5,3),
        (1,2,3,0),
        (0.7235, -1, -5, 0.3),
        (0.25, 1./3, 2, 0.999),
        (0.25, 1./3, 2, -1),
        (2,3,5,0.99),
        (3./2,-0.5,3,0.99),
        (2,2.5,-3.25,0.999),
        (-8, 18.016500331508873, 10.805295997850628, 0.90875647507000001),
        (-10,900,-10.5,0.99),
        (-10,900,10.5,0.99),
        (-1,2,1,1.0),
        (-1,2,1,-1.0),
        (-3,13,5,1.0),
        (-3,13,5,-1.0),
        (0.5, 1 - 270.5, 1.5, 0.999**2),  # from issue 1561
    ]
    dataset = [p + (float(mpmath.hyp2f1(*p)),) for p in pts]
    dataset = np.array(dataset, dtype=np.float_)

    olderr = np.seterr(invalid='ignore')
    try:
        FuncData(sc.hyp2f1, dataset, (0,1,2,3), 4, rtol=1e-10).check()
    finally:
        np.seterr(**olderr) 

def test_hyp2f1_real_some_points():
    pts = [
        (1, 2, 3, 0),
        (1./3, 2./3, 5./6, 27./32),
        (1./4, 1./2, 3./4, 80./81),
        (2,-2, -3, 3),
        (2, -3, -2, 3),
        (2, -1.5, -1.5, 3),
        (1, 2, 3, 0),
        (0.7235, -1, -5, 0.3),
        (0.25, 1./3, 2, 0.999),
        (0.25, 1./3, 2, -1),
        (2, 3, 5, 0.99),
        (3./2, -0.5, 3, 0.99),
        (2, 2.5, -3.25, 0.999),
        (-8, 18.016500331508873, 10.805295997850628, 0.90875647507000001),
        (-10, 900, -10.5, 0.99),
        (-10, 900, 10.5, 0.99),
        (-1, 2, 1, 1.0),
        (-1, 2, 1, -1.0),
        (-3, 13, 5, 1.0),
        (-3, 13, 5, -1.0),
        (0.5, 1 - 270.5, 1.5, 0.999**2),  # from issue 1561
    ]
    dataset = [p + (float(mpmath.hyp2f1(*p)),) for p in pts]
    dataset = np.array(dataset, dtype=np.float_)

    olderr = np.seterr(invalid='ignore')
    try:
        FuncData(sc.hyp2f1, dataset, (0,1,2,3), 4, rtol=1e-10).check()
    finally:
        np.seterr(**olderr) 

def test_hyp2f1(self):
        assert_mpmath_equal(sc.hyp2f1,
                            _exception_to_nan(lambda a, b, c, x: mpmath.hyp2f1(a, b, c, x, **HYPERKW)),
                            [Arg(), Arg(), Arg(), Arg()])