Python scipy.integrate 模块，dblquad() 实例源码

def _B_1_function(self, z, B_0=None):
        """
        calculate B_1(z) function defined in:

        Gould A. 1994 ApJ 421L, 71 "Proper motions of MACHOs
        http://adsabs.harvard.edu/abs/1994ApJ...421L..71G

        Yoo J. et al. 2004 ApJ 603, 139 "OGLE-2003-BLG-262: Finite-Source
        Effects from a Point-Mass Lens"
        http://adsabs.harvard.edu/abs/2004ApJ...603..139Y

        """
        if B_0 is None:
            B_0 = self._B_0_function(z)

        function = (lambda r, theta: r * np.sqrt(1.-r**2) /
            self.parameters.rho / np.sqrt(r**2+zz**2-2.*r*zz*np.cos(theta)))
        lim_0 = lambda x: 0
        lim_1 = lambda x: 1
        W_1 = 0. * z
        for (i, zz) in enumerate(z):
            W_1[i] = integrate.dblquad(function, 0., 2.*np.pi, lim_0, lim_1)[0]

        W_1 /= np.pi
        return B_0 - 1.5 * z * self.parameters.rho * W_1

def z(self):
        """
        Calculate the partion function .
        """
        from scipy.integrate import dblquad
        from numpy import pi

        return dblquad(self.prob, 0., 2 * pi, lambda x: 0., lambda x: 2 * pi)

def test_finite_line_source(self, rel_tol=1.0e-6):
        """ Tests the value of the FLS solution.
        """
        from pygfunction.boreholes import Borehole
        from pygfunction.heat_transfer import finite_line_source
        # Evaluate the double integral
        reference = dblquad(fls_double,
                            self.D1, self.D1+self.H1,
                            lambda x: self.D2, lambda x: self.D2+self.H2,
                            args=(self.t, self.dis, self.alpha))[0]/self.H2
        # Evaluate using heat_transfer.finite_line_source
        borehole1 = Borehole(self.H1, self.D1, 0.05, 0., 0.)
        borehole2 = Borehole(self.H2, self.D2, 0.05, self.dis, 0.)
        calculated = finite_line_source(self.t, self.alpha,
                                        borehole1, borehole2)
        self.assertAlmostEqual(calculated, reference,
                               delta=rel_tol*reference,
                               msg='Incorrect value of finite line source '
                                   'solution.')

def test_finite_line_source_real_part(self, rel_tol=1.0e-6):
        """ Tests the value of the real part of the FLS solution.
        """
        from pygfunction.boreholes import Borehole
        from pygfunction.heat_transfer import finite_line_source
        # Evaluate the double integral
        reference = dblquad(fls_double,
                            self.D1, self.D1+self.H1,
                            lambda x: self.D2, lambda x: self.D2+self.H2,
                            args=(self.t,
                                  self.dis,
                                  self.alpha,
                                  True,
                                  False))[0]/self.H2
        # Evaluate using heat_transfer.finite_line_source
        borehole1 = Borehole(self.H1, self.D1, 0.05, 0., 0.)
        borehole2 = Borehole(self.H2, self.D2, 0.05, self.dis, 0.)
        calculated = finite_line_source(self.t, self.alpha,
                                        borehole1, borehole2,
                                        reaSource=True, imgSource=False)
        self.assertAlmostEqual(calculated, reference,
                               delta=rel_tol*reference,
                               msg='Incorrect value of the real part of the '
                                   'finite line source solution.')

def test_finite_line_source_image_part(self, rel_tol=1.0e-6):
        """ Tests the value of the image part of the FLS solution.
        """
        from pygfunction.boreholes import Borehole
        from pygfunction.heat_transfer import finite_line_source
        # Evaluate the double integral
        reference = dblquad(fls_double,
                            self.D1, self.D1+self.H1,
                            lambda x: self.D2, lambda x: self.D2+self.H2,
                            args=(self.t,
                                  self.dis,
                                  self.alpha,
                                  False,
                                  True))[0]/self.H2
        # Evaluate using heat_transfer.finite_line_source
        borehole1 = Borehole(self.H1, self.D1, 0.05, 0., 0.)
        borehole2 = Borehole(self.H2, self.D2, 0.05, self.dis, 0.)
        calculated = finite_line_source(self.t, self.alpha,
                                        borehole1, borehole2,
                                        reaSource=False, imgSource=True)
        self.assertAlmostEqual(calculated, reference,
                               delta=np.abs(rel_tol*reference),
                               msg='Incorrect value of the image part of the '
                                   'finite line source solution.')

def test_price1(F1, F2, dv1, dv2, rho, K1, K2, T):
    v1 = dv1 * np.sqrt(T)
    v2 = dv2 * np.sqrt(T)
    var = scipy.stats.multivariate_normal(mean=[F1,F2], cov=[[v1*v2,v1*v2*rho],[v1*v2*rho,v2*v2]])
    def int_func1(x, y):
        return var.pdf([x,y])*(y-K2)

    def int_func2(x, y):
        return var.pdf([x,y])*(x-K1)

    res1 = dblquad(int_func1, K2, np.inf, lambda x: x-K2+K1, lambda x: np.inf)
    res2 = dblquad(int_func2, K2, np.inf, lambda x: K1, lambda x: x-K2+K1)
    return res1[0] + res2[0]

def mutual_information(self):
        def mi_integrand(x, y):
            return np.exp(self.logpdf_xy(x,y)) * \
                (self.logpdf_xy(x,y) - self.logpdf_x(x) - self.logpdf_y(y))
        return dblquad(
            lambda y, x: mi_integrand(x,y), self.D[0], self.D[1],
            self._lower_y, self._upper_y)

def _sanity_test(self):
        # Marginal of x integrates to one.
        print quad(lambda x: np.exp(self.logpdf_x(x)), self.D[0], self.D[1])

        # Marginal of y integrates to one.
        print quad(lambda y: np.exp(self.logpdf_y(y)), -1 ,1)

        # Joint of x,y integrates to one; quadrature will fail for small noise.
        print dblquad(
            lambda y,x: np.exp(self.logpdf_xy(x,y)), self.D[0], self.D[1],
            lambda x: -1, lambda x: 1)