Python integrate.quad函数代码示例

def Lookback(M1,L1,K1,function,p):
	if p == 1:
		file = open(function+'_'+str(M1)+'_'+str(L1)+'.txt','wt')
	LB=[]
	LB2=[]
	x=[]
	args = (M1,L1,K1)
	for z in drange(0,5,0.1):
	        x.append(z)
	        result, err = integrate.quad(E_z,0.0,z,args) # Lookback
	        result2, err = integrate.quad(E_z2,z,inf,args) # Age of Universe
	        LB.append(result)
	        LB2.append(result2)
		if p ==1:
	        	file.writelines(str(z)+" "+str(result)+" "+str(result2)+"\n")
	
	if p ==1:
		file.close
	plot(x,LB)
	plot(x,LB2)
	print LB[-1], LB2[-1]
	ylabel('Lookback Time $T_L/T_H$')
	xlabel('Redshift $z$')
	if p!= 1:
		show()

def foc(gov_policies, psi, sig, start=0, end=10):
    # Initialize local variables for government policies for better
    # readability.
    tax = gov_policies[0]
    trans = gov_policies[1]
    result = []

    # Compute different parts of first FOC (respect tax rate) and combine
    part_a = integrate.quad(
        lambda w: dell_u_tax(w, cons(w, tax, psi, trans),
                             hours(w, tax, psi), psi, tax) * lognorm.pdf(
                                 w, s=sig, scale=np.exp(-sig**2 / 2)),
        0, 10                   # set integration borders
    )[0]
    part_b = integrate.quad(
        lambda w: w * hours(w, tax, psi), start, end
    )[0]

    # Compute first part of the second FOC (respect transfers)
    part_c = integrate.quad(
        lambda w: lognorm.pdf(w, s=sig, scale=np.exp(-sig**2 / 2)) *
        dell_u_trans(cons(w, tax, psi, trans), hours(w, tax, psi), psi),
        start, end
    )[0]

    # Store first foc in results vector
    result.append(part_a + part_c * part_b)

    # Compute budget constraint
    bud_const = trans - integrate.quad(
        lambda w: tax * w * hours(w, tax, psi), start, end
    )[0]
    result.append(bud_const)

    return result

def _calculate_levy(x, alpha, beta, cdf=False):
    """ Calculation of Levy stable distribution via numerical integration.
        This is used in the creation of the lookup table. """
    # "0" parameterization as per http://academic2.american.edu/~jpnolan/stable/stable.html
    # Note: fails for alpha=1.0
    #       (so make sure alpha=1.0 isn't exactly on the interpolation grid)
    from scipy import integrate

    C = beta * N.tan(N.pi*0.5*alpha)
    
    def func_cos(u):
        ua = u ** alpha
        if ua > 700.0: return 0.0
        return N.exp(-ua)*N.cos(C*ua-C*u)

    def func_sin(u):
        ua = u ** alpha
        if ua > 700.0: return 0.0
        return N.exp(-ua)*N.sin(C*ua-C*u)

    if cdf:
        # Cumulative density function
        return (integrate.quad(lambda u: u and func_cos(u)/u or 0.0, 0.0, integrate.Inf, weight="sin", wvar=x, limlst=1000)[0]
               + integrate.quad(lambda u: u and func_sin(u)/u or 0.0, 0.0, integrate.Inf, weight="cos", wvar=x, limlst=1000)[0] 
               ) / N.pi + 0.5
    else:
        # Probability density function
        return ( integrate.quad(func_cos, 0.0, integrate.Inf, weight="cos", wvar=x, limlst=1000)[0]
               - integrate.quad(func_sin, 0.0, integrate.Inf, weight="sin", wvar=x, limlst=1000)[0] 
               ) / N.pi

    def radEinstein(self, zs):
        '''
        SDSSObject.radEinstein(z1, z2)
        ==============================
        Estimated Einstein Radius from Single Isothermal Sphere (SIS) model

        Parameters:
            zs: source redshift
        Returns:
            None
        '''
        # Necessary parameter
        H0 = 72e3
        c = 299792458.0
        OmegaM = 0.258
        OmegaL = 0.742
        vdisp = self.vdisp * 1000.0

        # Function needed for numerical computation of angular diameter distance
        def x12(z):
            return 1.0 / np.sqrt((1.0 - OmegaM - OmegaL) * (1.0 + z) *
                                 (1.0 + z) + OmegaM * (1.0 + z) ** 3.0 + OmegaL)
        # compute ang. diam. distances
        Ds = ((c / H0) * quad(x12, 0.0, zs)[0]) / (1.0 + zs)
        Dls = ((c / H0) * quad(x12, self.z, zs)[0]) / (1.0 + zs)
        # return value in arcsec
        coeff = 3600.0 * (180.0 / np.pi)
        return coeff * 4.0 * np.pi * vdisp * vdisp * Dls / (c * c * Ds)

def phi(d, dp):
    """vdW-DF kernel."""
    from scipy.integrate import quad

    kwargs = dict(epsabs=1.0e-6, epsrel=1.0e-6, limit=400)
    cut = 35
    return quad(lambda y: quad(f, 0, cut, (y, d, dp), **kwargs)[0], 0, cut, **kwargs)[0]

 def getConsumerWelfare(self, equilibriumValue):
     """
     Calculate total consumer welfare
         W = \int_{A} u(t,A) dt +\int_{B} u(t,B) dt.
     Can break this up as:
     W = Network Effect Benefits + Heterogeneity Stuff
     W = t^{e} * \nu(t^{e}, A) + (1-t^{e}) * \nu(1-t^{e}) + Heterogeneity
     NB: we are assuming demand has been normalized to 1
     
     @param equilibriumValue: the equilibrium to use when calculating social welfare
     @return: total consumer welfare
     """
     networkAPart = equilibriumValue * self._nu_A(equilibriumValue)
     # Remember _nu_B defined a function of network A size ...
     networkBPart = (1.0 - equilibriumValue) * \
         self._nu_B(equilibriumValue)
     
     hetAPart = integrate.quad(self._het_A, 0, equilibriumValue)[0]
     hetBPart = integrate.quad(self._het_B, equilibriumValue, 1)[0]
     
     priceAPart = equilibriumValue * self._priceFunction_A(self._price_A)
     priceBPart = equilibriumValue * self._priceFunction_B(self._price_B)
     
     # logger.debug('Welfare: (netA, netB, hetA, hetB): ' + \
     #    str((networkAPart, networkBPart, hetAPart, hetBPart)))
     totalWelfare = (networkAPart + networkBPart +
                     hetAPart + hetBPart + priceAPart + priceBPart)
     return totalWelfare

def complex_quad(func, a, b, **kw_args):
    func_re = lambda t: np.real(func(t))
    func_im = lambda t: np.imag(func(t))
    I_re = quad(func_re, a, b, **kw_args)[0]
    I_im = quad(func_im, a, b, **kw_args)[0]
    
    return I_re + 1j*I_im

    def test_b_less_than_a_3(self):
        def f(x):
            return 1.0

        val_1, err_1 = quad(f, 0, 1, weight='alg', wvar=(0, 0))
        val_2, err_2 = quad(f, 1, 0, weight='alg', wvar=(0, 0))
        assert_allclose(val_1, -val_2, atol=max(err_1, err_2))

    def test_ctypes_variants(self):
        sin_0 = get_clib_test_routine('_sin_0', ctypes.c_double,
                                      ctypes.c_double, ctypes.c_void_p)

        sin_1 = get_clib_test_routine('_sin_1', ctypes.c_double,
                                      ctypes.c_int, ctypes.POINTER(ctypes.c_double),
                                      ctypes.c_void_p)

        sin_2 = get_clib_test_routine('_sin_2', ctypes.c_double,
                                      ctypes.c_double)

        sin_3 = get_clib_test_routine('_sin_3', ctypes.c_double,
                                      ctypes.c_int, ctypes.POINTER(ctypes.c_double))

        sin_4 = get_clib_test_routine('_sin_3', ctypes.c_double,
                                      ctypes.c_int, ctypes.c_double)

        all_sigs = [sin_0, sin_1, sin_2, sin_3, sin_4]
        legacy_sigs = [sin_2, sin_4]
        legacy_only_sigs = [sin_4]

        # LowLevelCallables work for new signatures
        for j, func in enumerate(all_sigs):
            callback = LowLevelCallable(func)
            if func in legacy_only_sigs:
                assert_raises(ValueError, quad, callback, 0, pi)
            else:
                assert_allclose(quad(callback, 0, pi)[0], 2.0)

        # Plain ctypes items work only for legacy signatures
        for j, func in enumerate(legacy_sigs):
            if func in legacy_sigs:
                assert_allclose(quad(func, 0, pi)[0], 2.0)
            else:
                assert_raises(ValueError, quad, func, 0, pi)

    def test_b_less_than_a(self):
        def f(x, p, q):
            return p * np.exp(-q*x)

        val_1, err_1 = quad(f, 0, np.inf, args=(2, 3))
        val_2, err_2 = quad(f, np.inf, 0, args=(2, 3))
        assert_allclose(val_1, -val_2, atol=max(err_1, err_2))

    def test_b_less_than_a_2(self):
        def f(x, s):
            return np.exp(-x**2 / 2 / s) / np.sqrt(2.*s)

        val_1, err_1 = quad(f, -np.inf, np.inf, args=(2,))
        val_2, err_2 = quad(f, np.inf, -np.inf, args=(2,))
        assert_allclose(val_1, -val_2, atol=max(err_1, err_2))

  def get_FO(self,x,y,z,Q2,qT2,muR2,muF2,charge,ps='dxdQ2dzdqT2',method='gauss'):
    D=self.D
    D['A1']=1+(2/y-1)**2
    D['A2']=-2
    w2=qT2/Q2*x*z
    w=w2**0.5
    xia_=lambda xib: w2/(xib-z)+x
    xib_=lambda xia: w2/(xia-x)+z

    integrand_xia=lambda xia: self.get_M(xia,xib_(xia),x/xia,z/xib_(xia),Q2,muF2,qT2,charge)
    integrand_xib=lambda xib: self.get_M(xia_(xib),xib,x/xia_(xib),z/xib,Q2,muF2,qT2,charge)

    if method=='quad':
      FO = quad(integrand_xia,x+w,1)[0] + quad(integrand_xib,z+w,1)[0]

    elif method=='gauss':
      integrand_xia=np.vectorize(integrand_xia)
      integrand_xib=np.vectorize(integrand_xib)
      FO = fixed_quad(integrand_xia,x+w,1,n=40)[0] + fixed_quad(integrand_xib,z+w,1,n=40)[0]

    if ps=='dxdQ2dzdqT2':
      s=x*y*Q2
      prefactor = D['alphaEM']**2 * self.SC.get_alphaS(muR2) 
      prefactor/= 2*s**2*Q2*x**2
      prefactor*= D['GeV**-2 -> nb'] 
      return prefactor * FO
    else: 
      print 'ps not inplemented'
      return None

def nc_of_norm():
    f1 = lambda x: x**4
    f2 = lambda x: x**2-x+2
    
    rv = norm(loc = 2 , scale=20)
    print rv.mean()
    print rv.var()
    print rv.moment(1)
    print rv.moment(4)
    #     moments的参数可为m(均值),v(方差值),s(偏度),k(峰度),默认为mv
    print rv.stats(moments = 'mvsk')
# 3）scipy.stats.norm.expect(func,loc=0,scale=1)函数返回func函数相对于正态分布的期望值，其中函数f(x)相对于分布dist的期望值定义为E[x]=Integral(f(x) * dist.pdf(x))
    print(norm.expect(f1, loc=1, scale=2))
    print(norm.expect(f2, loc=2, scale=5))
    
    
    
# （2）lambda argument_list:expression用来表示匿名函数
# （3）numpy.exp(x)计算x的指数
# （4）numpy.inf表示正无穷大
# （5）scipy.integrate.quad(func,a,b)计算func函数从a至b的积分
# （6）scipy.stats.expon(scale)函数返回符合指数分布的参数为scale的随机变量rv
# （7）rv.moment(n,*args,*kwds)返回随机变量的n阶距
        #1st non-center moment of expon distribution whose lambda is 0.5
    E1 = lambda x: x*0.5*np.exp(-x/2)
    #2nd non-center moment of expon distribution whose lambda is 0.5
    E2 = lambda x: x**2*0.5*np.exp(-x/2)
    print(integrate.quad(E1, 0, np.inf))
    print(integrate.quad(E2, 0, np.inf))

 def computeColor(self, filtX, filtY, z):
     """Compute color (flux in filter X - filter Y) of SED at redshift z, return color in magnitudes
     
        @param filtX   lower wavelength filter
        @param filtY   higher wavelength filter
        @param z       redshift of SED
     """
     
     if filtX not in self.filterDict:
         emsg = "Filter " + filtX + " is not in the filter dictionary"
         raise LookupError(emsg)
     if filtY not in self.filterDict:
         emsg = "Filter " + filtY + " is not in the filter dictionary"
         raise LookupError(emsg)
     if filtX == filtY:
         raise ValueError("ERROR! cannot have color as difference between same filter")
     
     # define integral limits by filter extents
     aX, bX = self.filterDict[filtX].returnFilterRange()
     aY, bY = self.filterDict[filtY].returnFilterRange()
     
     int1 = integ.quad(self._integrand1, aX, bX, args=(z, filtX))[0]
     int2 = integ.quad(self._integrand1, aX, bX, args=(z, filtY))[0]
     
     # Do we need this zero-point term?
     int3 = integ.quad(self._integrand2, aX, bX, args=(filtX))[0]
     int4 = integ.quad(self._integrand2, aX, bX, args=(filtY))[0]
     zp = -2.5*math.log10(int4/int3)
     
     return -2.5*math.log10(int1/int2) + zp

 def funcf(E,verbose=False):
     epsilon = 0
     try:
         t = E.shape
         fans = []
         problems = []
         for i in range(len(E)):
             print i+1, ' of ', len(E)
             rapoval = rapo(E[i])
             prefactor = (1./(sqrt(8)*pi**2*(model.Mnorm + 10**Mencgood(log10(rapoval)))))
             temp = intg.quad(finterior,epsilon,1-epsilon,args=(E[i],rapoval),full_output = 1)
             t = temp[0]
             try:
                 if temp[3] != '':
                     if verbose == True:
                         print 'f, E = ',E[i],'message = ',temp[3]
                     problems.append(i)                  
             except IndexError:
                 pass
             fans.append((prefactor*t)[0])
         return array(fans),problems
     except AttributeError:
         rapoval = rapo(E)
         prefactor = (1./(sqrt(8)*pi**2*(model.Mnorm + 10**Mencgood(log10(rapoval)))))
         temp = intg.quad(finterior,epsilon,1-epsilon,args=(E,rapoval),full_output = 1)
         t = temp[0]
         problem = []
         try:
             if temp1[3] != '':
                 if verbose == True:
                     print 'f, E = ',E,'message = ',temp1[3]
                 problem = [E]
         except IndexError:
             pass
         return prefactor*t, problem