Python scipy.integrate.nquad() Examples

def getThresholdMoments(x,myP,whichModel):
    if whichModel==0: # Gaussian
        M1,err = nInt.quad(integrateGaussianMoment,-5,5,args=(x[0],myP,1)) 
        M2,err = nInt.quad(integrateGaussianMoment,-5,5,args=(x[0],myP,2)) 
    elif whichModel==1: # t
        lowerBound = np.maximum(x[1]-40,2)
        support = [[-7,7],[lowerBound,x[1]+40]]
        M1,err=nInt.nquad(thresholdMoment,support,args=(x[0],x[1],myP,whichModel,1))
        M2,err=nInt.nquad(thresholdMoment,support,args=(x[0],x[1],myP,whichModel,2))
    elif whichModel==2: # Variance-gamma
        invCdf = th.nvmPpf(myP,x[1],0)
        support = [[-7,7],[0,100]]
        M1,err=nInt.nquad(thresholdMoment,support,args=(x[0],x[1],myP,whichModel,1,invCdf))
        M2,err=nInt.nquad(thresholdMoment,support,args=(x[0],x[1],myP,whichModel,2,invCdf))
    elif whichModel==3: # Generalized hyperbolic        
        invCdf = th.nvmPpf(myP,x[1],1)
        support = [[-7,7],[0,100]]
        M1,err=nInt.nquad(thresholdMoment,support,args=(x[0],x[1],myP,whichModel,1,invCdf))
        M2,err=nInt.nquad(thresholdMoment,support,args=(x[0],x[1],myP,whichModel,2,invCdf))
    return M1,M2 

def test_square_aliased_ranges_and_opts(self):
        def f(y, x):
            return 1.0

        r = [-1, 1]
        opt = {}
        assert_quad(nquad(f, [r, r], opts=[opt, opt]), 4.0) 

def test(self, function):
        area, _ = nquad(function, self.space)
        passed = abs(1 - area) <= self.epsilon
        self.areas.append(area)
        self.results.append(passed)

        return passed 

def jointDefaultProbabilityNVM(p,q,myRho,myA,whichModel):    
    invCdf = nvmPpf(p,myA,whichModel)
    support = [[-8,8],[0,100]]
    pr,err=nInt.nquad(jointIntegrandNVM,support,args=(p,q,myRho,myA,invCdf,whichModel))
    return pr 

def jointDefaultProbabilityT(p,q,myRho,nu):
    lowerBound = np.maximum(nu-40,2)
    support = [[-10,10],[lowerBound,nu+40]]
    pr,err=nInt.nquad(jointIntegrandT,support,args=(p,q,myRho,nu))
    return pr 

def myESCYT(l,p,c,p1,p2,whichModel):
    lowerBound = np.maximum(p2-20,0.0001)
    support = [[-8,8],[lowerBound,p2+8]]
    myAlpha = myApproxT(l,p,c,p1,p2,whichModel,1)
    esC = np.zeros(len(p))
    for n in range(0,len(p)):
        esC[n],err = nInt.nquad(integrateAllT,support,args=(l,n,p,c,p1,p2,whichModel))
    return (1/myAlpha)*esC 

def myVaRCYT(l,p,c,p1,p2,whichModel):
    lowerBound = np.maximum(p2-20,0.0001)
    support = [[-8,8],[lowerBound,p2+8]]
    den = myApproxT(l,p,c,p1,p2,whichModel,0)
    num = np.zeros(len(p))
    for n in range(0,len(p)):
        num[n],err = nInt.nquad(varCNumeratorT,support,args=(l,n,p,c,p1,p2,whichModel))
    return c*np.divide(num,den) 

def myApproxT(l,p,c,p1,p2,whichModel,myDegree,constant=0,den=1):
    lowerBound = np.maximum(p2-20,0.0001)
    support = [[-8,8],[lowerBound,p2+8]]
    if myDegree==2:
        constant = l
        den,err = nInt.nquad(computeYIntegralT,support,args=(l,p,c,p1,p2,whichModel,1))        
    num,err = nInt.nquad(computeYIntegralT,support,args=(l,p,c,p1,p2,whichModel,myDegree))
    return constant + np.divide(num,den) 

def test_gauss2darea(self):
        p = numpy.copy(self.p2d)
        p[2] = self.sigma
        p[3] = (numpy.random.rand() + 0.1) * self.sigma
        area = xu.math.Gauss2dArea(*p)
        (numarea, err) = nquad(xu.math.Gauss2d, [[-numpy.inf, numpy.inf],
                                                 [-numpy.inf, numpy.inf]],
                               args=tuple(p))
        digits = int(numpy.abs(numpy.log10(err))) - 3
        self.assertTrue(digits >= 3)
        self.assertAlmostEqual(area, numarea, places=digits) 

def test_scipy():
    from scipy import integrate
    import numpy as np
    args = t, x = se.symbols('t, x')
    lmb = se.Lambdify(args, [se.exp(-x*t)/t**5], as_scipy=True)
    res = integrate.nquad(lmb, [[1, np.inf], [0, np.inf]])
    assert abs(res[0] - 0.2) < 1e-7 

def test_dict_as_opts(self):
        try:
            out = nquad(lambda x, y: x * y, [[0, 1], [0, 1]], opts={'epsrel': 0.0001})
        except(TypeError):
            assert False 

def test_matching_tplquad(self):
        def func3d(x0, x1, x2, c0, c1):
            return x0**2 + c0 * x1**3 - x0 * x1 + 1 + c1 * np.sin(x2)

        res = tplquad(func3d, -1, 2, lambda x: -2, lambda x: 2,
                      lambda x, y: -np.pi, lambda x, y: np.pi,
                      args=(2, 3))
        res2 = nquad(func3d, [[-np.pi, np.pi], [-2, 2], (-1, 2)], args=(2, 3))
        assert_almost_equal(res, res2) 

def test_matching_dblquad(self):
        def func2d(x0, x1):
            return x0**2 + x1**3 - x0 * x1 + 1

        res, reserr = dblquad(func2d, -2, 2, lambda x: -3, lambda x: 3)
        res2, reserr2 = nquad(func2d, [[-3, 3], (-2, 2)])
        assert_almost_equal(res, res2)
        assert_almost_equal(reserr, reserr2) 

def test_matching_quad(self):
        def func(x):
            return x**2 + 1

        res, reserr = quad(func, 0, 4)
        res2, reserr2 = nquad(func, ranges=[[0, 4]])
        assert_almost_equal(res, res2)
        assert_almost_equal(reserr, reserr2) 

def test_square_aliased_fn_ranges_and_opts(self):
        def f(y, x):
            return 1.0

        def fn_range(*args):
            return (-1, 1)

        def fn_opt(*args):
            return {}

        ranges = [fn_range, fn_range]
        opts = [fn_opt, fn_opt]
        assert_quad(nquad(f, ranges, opts=opts), 4.0) 

def test_square_separate_ranges_and_opts(self):
        def f(y, x):
            return 1.0

        assert_quad(nquad(f, [[-1, 1], [-1, 1]], opts=[{}, {}]), 4.0) 

def test_variable_limits(self):
        scale = .1

        def func2(x0, x1, x2, x3, t0, t1):
            val = (x0*x1*x3**2 + np.sin(x2) + 1 +
                   (1 if x0 + t1*x1 - t0 > 0 else 0))
            return val

        def lim0(x1, x2, x3, t0, t1):
            return [scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) - 1,
                    scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) + 1]

        def lim1(x2, x3, t0, t1):
            return [scale * (t0*x2 + t1*x3) - 1,
                    scale * (t0*x2 + t1*x3) + 1]

        def lim2(x3, t0, t1):
            return [scale * (x3 + t0**2*t1**3) - 1,
                    scale * (x3 + t0**2*t1**3) + 1]

        def lim3(t0, t1):
            return [scale * (t0 + t1) - 1, scale * (t0 + t1) + 1]

        def opts0(x1, x2, x3, t0, t1):
            return {'points': [t0 - t1*x1]}

        def opts1(x2, x3, t0, t1):
            return {}

        def opts2(x3, t0, t1):
            return {}

        def opts3(t0, t1):
            return {}

        res = nquad(func2, [lim0, lim1, lim2, lim3], args=(0, 0),
                    opts=[opts0, opts1, opts2, opts3])
        assert_quad(res, 25.066666666666663) 

def test_fixed_limits(self):
        def func1(x0, x1, x2, x3):
            val = (x0**2 + x1*x2 - x3**3 + np.sin(x0) +
                   (1 if (x0 - 0.2*x3 - 0.5 - 0.25*x1 > 0) else 0))
            return val

        def opts_basic(*args):
            return {'points': [0.2*args[2] + 0.5 + 0.25*args[0]]}

        res = nquad(func1, [[0, 1], [-1, 1], [.13, .8], [-.15, 1]],
                    opts=[opts_basic, {}, {}, {}], full_output=True)
        assert_quad(res[:-1], 1.5267454070738635)
        assert_(res[-1]['neval'] > 0 and res[-1]['neval'] < 4e5) 

def test_integrate_2d(self):
        np.random.seed(1234)
        c = np.random.rand(4, 5, 16, 17)
        x = np.linspace(0, 1, 16+1)**1
        y = np.linspace(0, 1, 17+1)**2

        # make continuously differentiable so that nquad() has an
        # easier time
        c = c.transpose(0, 2, 1, 3)
        cx = c.reshape(c.shape[0], c.shape[1], -1).copy()
        _ppoly.fix_continuity(cx, x, 2)
        c = cx.reshape(c.shape)
        c = c.transpose(0, 2, 1, 3)
        c = c.transpose(1, 3, 0, 2)
        cx = c.reshape(c.shape[0], c.shape[1], -1).copy()
        _ppoly.fix_continuity(cx, y, 2)
        c = cx.reshape(c.shape)
        c = c.transpose(2, 0, 3, 1).copy()

        # Check integration
        p = NdPPoly(c, (x, y))

        for ranges in [[(0, 1), (0, 1)],
                       [(0, 0.5), (0, 1)],
                       [(0, 1), (0, 0.5)],
                       [(0.3, 0.7), (0.6, 0.2)]]:

            ig = p.integrate(ranges)
            ig2, err2 = nquad(lambda x, y: p((x, y)), ranges,
                              opts=[dict(epsrel=1e-5, epsabs=1e-5)]*2)
            assert_allclose(ig, ig2, rtol=1e-5, atol=1e-5,
                            err_msg=repr(ranges)) 

def ise(estimateddensity, theoreticaldensity):
    def integrand(*x):
        value = np.array(x)
        return (estimateddensity(value) - theoreticaldensity(value))**2

    return integrate.nquad(integrand, [[-5., 5.], [-5., 5.]]) 

def test_matching_tplquad(self):
        def func3d(x0, x1, x2, c0, c1):
            return x0**2 + c0 * x1**3 - x0 * x1 + 1 + c1 * np.sin(x2)


        res = tplquad(func3d, -1, 2, lambda x:-2, lambda x:2,
                      lambda x, y : -np.pi, lambda x, y : np.pi,
                      args=(2, 3))
        res2 = nquad(func3d, [[-np.pi, np.pi], [-2, 2], (-1, 2)], args=(2, 3))
        assert_almost_equal(res, res2) 

def test_matching_dblquad(self):
        def func2d(x0, x1):
            return x0**2 + x1**3 - x0 * x1 + 1

        res, reserr = dblquad(func2d, -2, 2, lambda x:-3, lambda x:3)
        res2, reserr2 = nquad(func2d, [[-3, 3], (-2, 2)])
        assert_almost_equal(res, res2)
        assert_almost_equal(reserr, reserr2) 

def test_matching_quad(self):
        def func(x):
            return x**2 + 1

        res, reserr = quad(func, 0, 4)
        res2, reserr2 = nquad(func, ranges=[[0, 4]])
        assert_almost_equal(res, res2)
        assert_almost_equal(reserr, reserr2) 

def test_square_separate_fn_ranges_and_opts(self):
        def f(y, x):
            return 1.0
        def fn_range0(*args):
            return (-1, 1)
        def fn_range1(*args):
            return (-1, 1)
        def fn_opt0(*args):
            return {}
        def fn_opt1(*args):
            return {}

        ranges = [fn_range0, fn_range1]
        opts = [fn_opt0, fn_opt1]
        assert_quad(nquad(f, ranges, opts=opts), 4.0) 

def test_square_aliased_ranges_and_opts(self):
        def f(y, x):
            return 1.0

        r = [-1, 1]
        opt = {}
        assert_quad(nquad(f, [r, r], opts=[opt, opt]), 4.0) 

def test_square_separate_ranges_and_opts(self):
        def f(y, x):
            return 1.0

        assert_quad(nquad(f, [[-1, 1], [-1, 1]], opts=[{}, {}]), 4.0) 

def test_variable_limits(self):
        scale = .1
        def func2(x0, x1, x2, x3, t0, t1):
            return x0*x1*x3**2 + np.sin(x2) + 1 + (1 if x0 + t1*x1 - t0 > 0 else 0)
        def lim0(x1, x2, x3, t0, t1):
            return [scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) - 1,
                    scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) + 1]
        def lim1(x2, x3, t0, t1):
            return [scale * (t0*x2 + t1*x3) - 1,
                    scale * (t0*x2 + t1*x3) + 1]
        def lim2(x3, t0, t1):
            return [scale * (x3 + t0**2*t1**3) - 1,
                    scale * (x3 + t0**2*t1**3) + 1]
        def lim3(t0, t1):
            return [scale * (t0 + t1) - 1, scale * (t0 + t1) + 1]
        def opts0(x1, x2, x3, t0, t1):
            return {'points':[t0 - t1*x1]}
        def opts1(x2, x3, t0, t1):
            return {}
        def opts2(x3, t0, t1):
            return {}
        def opts3(t0, t1):
            return {}

        res = nquad(func2, [lim0, lim1, lim2, lim3], args=(0,0),
                    opts=[opts0, opts1, opts2, opts3])
        assert_quad(res, 25.066666666666663) 

def test_fixed_limits(self):
        def func1(x0, x1, x2, x3):
            return x0**2 + x1*x2 - x3**3 + np.sin(x0) + (1 if
                        (x0 - 0.2*x3 - 0.5 - 0.25*x1 > 0) else 0)

        def opts_basic(*args):
            return {'points' : [0.2*args[2] + 0.5 + 0.25*args[0]]}

        res = nquad(func1, [[0,1], [-1,1], [.13,.8], [-.15,1]],
                    opts=[opts_basic, {}, {}, {}])
        assert_quad(res, 1.5267454070738635) 

def compute_expectation(model, fun):
    # Computes the expectation of fun(th) wrt the posterior distribution,
    # given by exp(model.log_p_marg(th))
    th_shape = model.th_shape
    L = np.prod(th_shape)
    lims = [[-np.inf, np.inf]] * L
    def integrand(*args):
        th = np.array(args).reshape(th_shape)
        return np.exp(model.log_p_marg(th)) * fun(th)

    return integrate.nquad(integrand, lims)[0] 

def dblquad(func, a, b, gfun, hfun, args=(), epsabs=1.49e-8, epsrel=1.49e-8):
    """
    Compute a double integral.

    Return the double (definite) integral of ``func(y, x)`` from ``x = a..b``
    and ``y = gfun(x)..hfun(x)``.

    Parameters
    ----------
    func : callable
        A Python function or method of at least two variables: y must be the
        first argument and x the second argument.
    a, b : float
        The limits of integration in x: `a` < `b`
    gfun : callable
        The lower boundary curve in y which is a function taking a single
        floating point argument (x) and returning a floating point result: a
        lambda function can be useful here.
    hfun : callable
        The upper boundary curve in y (same requirements as `gfun`).
    args : sequence, optional
        Extra arguments to pass to `func`.
    epsabs : float, optional
        Absolute tolerance passed directly to the inner 1-D quadrature
        integration. Default is 1.49e-8.
    epsrel : float, optional
        Relative tolerance of the inner 1-D integrals. Default is 1.49e-8.

    Returns
    -------
    y : float
        The resultant integral.
    abserr : float
        An estimate of the error.

    See also
    --------
    quad : single integral
    tplquad : triple integral
    nquad : N-dimensional integrals
    fixed_quad : fixed-order Gaussian quadrature
    quadrature : adaptive Gaussian quadrature
    odeint : ODE integrator
    ode : ODE integrator
    simps : integrator for sampled data
    romb : integrator for sampled data
    scipy.special : for coefficients and roots of orthogonal polynomials

    """
    def temp_ranges(*args):
        return [gfun(args[0]), hfun(args[0])]
    return nquad(func, [temp_ranges, [a, b]], args=args, 
            opts={"epsabs": epsabs, "epsrel": epsrel})