Python scipy.sqrt方法代码示例

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import sqrt [as 别名]
def _fug(self, Z, xi):
        rho=self.P.atm/Z/R_atml/self.T
        tita=[]
        for i in range(len(self.componente)):
            suma=0
            for j in range(len(self.componente)):
                suma+=xi[j]*(-self.Ao**0.5*self.Aoi[i]**0.5*(1-self.kij[i][j]) \
                                    -  self.Co**0.5*self.Coi[i]**0.5*(1-self.kij[i][j])**3/self.T**2 \
                                    + self.Do**0.5*self.Doi[i]**0.5*(1-self.kij[i][j])**4/self.T**3 \
                                    -  self.Eo**0.5*self.Eoi[i]**0.5*(1-self.kij[i][j])**5/self.T**4)
#                print suma
            lo=R_atml*self.T*log(rho*R_atml*self.T*xi[i]) \
                    + rho*(self.Bo+self.Boi[i])*R_atml*self.T \
                    + 2*rho*suma \
                    + rho**2/2*(3*(self.b**2*self.bi[i])**(1./3)*R_atml*self.T-3*(self.a**2*self.ai[i])**(1./3)-3*(self.d**2*self.di[i])**(1./3)/self.T) \
                    + self.alfa*rho**5/5*(3*(self.a**2*self.ai[i])**(1./3)+3*(self.d**2*self.di[i])**(1./3)/self.T) \
                    + 3*rho**5/5*(self.a+self.d/self.T)*(self.alfa**2*self.alfai[i])**(1./3) \
                    + 3*(self.c**2*self.ci[i])**(1./3)*rho**2/self.T**2*((1-exp(-self.gamma*rho**2))/self.gamma/rho**2-exp(-self.gamma*rho**2)/2) \
                    - (2*self.c*sqrt(self.gammai[i]/self.gamma)**0.5/self.gamma/self.T**2)*((1-exp(-self.gamma*rho**2))*(1+self.gamma*rho**2+self.gamma**2*rho**4/2))
            tita.append(exp(lo/R_atml/self.T))
        return tita 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import sqrt [as 别名]
def up_and_out_call(s0,x,T,r,sigma,n_simulation,barrier):
    n_steps=100. 
    dt=T/n_steps 
    total=0 
    for j in sp.arange(0, n_simulation): 
        sT=s0 
        out=False
        for i in range(0,int(n_steps)): 
            e=sp.random.normal() 
            sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt)) 
            if sT>barrier: 
               out=True 
        if out==False: 
            total+=bsCall(s0,x,T,r,sigma) 
    return total/n_simulation 
# 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import sqrt [as 别名]
def implied_vol_call_min(S,X,T,r,c): 
    from scipy import log,exp,sqrt,stats 
    implied_vol=1.0
    min_value=1000
    for i in range(10000): 
        sigma=0.0001*(i+1)
        d1=(log(S/X)+(r+sigma*sigma/2.)*T)/(sigma*sqrt(T)) 
        d2 = d1-sigma*sqrt(T)
        c2=S*stats.norm.cdf(d1)-X*exp(-r*T)*stats.norm.cdf(d2) 
        abs_diff=abs(c2-c)
        if abs_diff<min_value: 
            min_value=abs_diff 
            implied_vol=sigma 
            k=i
    return implied_vol


# Step 3: get call option data 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import sqrt [as 别名]
def lookback_min_price_as_strike(s,T,r,sigma,n_simulation): 
    n_steps=100
    dt=T/n_steps
    total=0
    for j in range(n_simulation): 
        min_price=100000.	# a very big number 
        sT=s
        for i in range(int(n_steps)): 
            e=sp.random.normal()
            sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt)) 
            if sT<min_price:
                min_price=sT
                #print 'j=',j,'i=',i,'total=',total 
                total+=p4f.bs_call(s,min_price,T,r,sigma)
    return total/n_simulation

# 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import sqrt [as 别名]
def up_and_out_call(s0,x,T,r,sigma,n_simulation,barrier):
    n_steps=100.
    dt=T/n_steps 
    total=0
    for j in range(0, n_simulation): 
        sT=s0
        out=False
        for i in range(0,int(n_steps)): 
            e=sp.random.normal()
            sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt)) 
            if sT>barrier:
                out=True 
    if out==False:
       total+=p4f.bs_call(s0,x,T,r,sigma) 
    return total/n_simulation
# 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import sqrt [as 别名]
def KMV_f(E,D,T,r,sigmaE):
    n=10000
    m=2000
    diffOld=1e6     # a very big number
    for i in sp.arange(1,10):
        for j in sp.arange(1,m):
            A=E+D/2+i*D/n
            sigmaA=0.05+j*(1.0-0.001)/m
            d1 = (log(A/D)+(r+sigmaA*sigmaA/2.)*T)/(sigmaA*sqrt(T))
            d2 = d1-sigmaA*sqrt(T)
            diff4E= (A*N(d1)-D*exp(-r*T)*N(d2)-E)/A  # scale by assets
            diff4A= A/E*N(d1)*sigmaA-sigmaE          # a small number already
            diffNew=abs(diff4E)+abs(diff4A)
            if diffNew<diffOld:
               diffOld=diffNew
               output=(round(A,2),round(sigmaA,4),round(diffNew,5))
    return output
# 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import sqrt [as 别名]
def _column_covariances(X, uniformity_thresh):
    """."""

    Xvert = _high_frequency_vert(X, sigma=4.0)
    Xvertp = _high_frequency_vert(X, sigma=3.0)
    models = []
    use_C = []
    for i in range(X.shape[2]):
        xsub = Xvert[:, :, i]
        xsubp = Xvertp[:, :, i]
        mu = xsub.mean(axis=0)
        dists = s.sqrt(pow((xsub - mu), 2).sum(axis=1))
        distsp = s.sqrt(pow((xsubp - mu), 2).sum(axis=1))
        thresh = _percentile(dists, 95.0)
        uthresh = dists * uniformity_thresh
        #use       = s.logical_and(dists<thresh, abs(dists-distsp) < uthresh)
        use = dists < thresh
        C = s.cov(xsub[use, :], rowvar=False)
        [U, V, D] = svd(C)
        V[V < 1e-8] = 1e-8
        C = U.dot(s.diagflat(V)).dot(D)
        models.append(C)
        use_C.append(use)
    return s.array(models), Xvert, Xvertp, s.array(use_C).T 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import sqrt [as 别名]
def generate_data(N, S, L):

        # generate genetics
        G = 1.0 * (sp.rand(N, S) < 0.2)
        G -= G.mean(0)
        G /= G.std(0) * sp.sqrt(G.shape[1])

        # generate latent phenotypes
        Zg = sp.dot(G, sp.randn(G.shape[1], L))
        Zn = sp.randn(N, L)

        # generate variance exapleind
        vg = sp.linspace(0.8, 0, L)

        # rescale and sum
        Zg *= sp.sqrt(vg / Zg.var(0))
        Zn *= sp.sqrt((1 - vg) / Zn.var(0))
        Z = Zg + Zn

        return Z, G 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import sqrt [as 别名]
def ldpred_fast(betas_blup, h2=None, Ns=None, p_c=None):
    """
    LDpred-fast
    
    Joel Mefford's sparsified BLUP shrink. 
    """
    M = len(betas_blup)
    h2mn = h2+M/Ns
    h2mpn = h2+(M*p_c)/Ns
    C1 = h2mn/h2mpn
    C2 = h2/h2mn
    V_nc = (C2/Ns)*(1-C2*h2)
    V_c = C2*(h2/(M*p_c))+V_nc
    P_beta_ncaus = stats.norm.pdf(betas_blup,scale=sp.sqrt(V_nc))
    P_beta_caus = stats.norm.pdf(betas_blup,scale=sp.sqrt(V_c))
    P_caus_beta =  P_beta_caus *p_c/(P_beta_caus *p_c + P_beta_ncaus *(1-p_c))
    beta_jms = betas_blup * C1 * P_caus_beta 
    return beta_jms 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import sqrt [as 别名]
def get_beta(pval_read, raw_beta, beta_read, line_dict, header_dict, se, 
             z_from_se, N, eff_type, se_inferred_zscores):
    if eff_type=='BLUP':
        return raw_beta
    if z_from_se:
        assert se in header_dict, "SE was not specified in summary statistics provided, which is necessary when the 'z_from_se' flag is used."
        se_read = float(line_dict[header_dict[se]])
        if se_read==0 or not isfinite(se_read):
            return None
        se_inferred_zscores += 1
        return get_beta_from_se(beta_read, se_read, eff_type, raw_beta, N)
    else:
        if pval_read==0 or not isfinite(stats.norm.ppf(pval_read)):
            #Attempt to Parse SEs to infer Z-score
            if not se in header_dict:
                return None
            else:
                se_read = float(line_dict[header_dict[se]])
                if se_read==0 or not isfinite(se_read):
                    return None

                se_inferred_zscores += 1
                return get_beta_from_se(beta_read, se_read, eff_type, raw_beta, N)
        else:               
            return sp.sign(raw_beta) * stats.norm.ppf(pval_read / 2.0)/ sp.sqrt(N) 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import sqrt [as 别名]
def compute_residuals(self):
        """Compute residuals and stopping thresholds."""

        if self.opt['AutoRho', 'StdResiduals']:
            r = linalg.norm(self.rsdl_r(self.AXnr, self.Y))
            s = linalg.norm(self.rsdl_s(self.Yprev, self.Y))
            epri = scipy.sqrt(self.Nc)*self.opt['AbsStopTol'] + \
                self.rsdl_rn(self.AXnr, self.Y)*self.opt['RelStopTol']
            edua = scipy.sqrt(self.Nx)*self.opt['AbsStopTol'] + \
                self.rsdl_sn(self.U)*self.opt['RelStopTol']
        else:
            rn = self.rsdl_rn(self.AXnr, self.Y)
            if rn == 0.0:
                rn = 1.0
            sn = self.rsdl_sn(self.U)
            if sn == 0.0:
                sn = 1.0
            r = linalg.norm(self.rsdl_r(self.AXnr, self.Y)) / rn
            s = linalg.norm(self.rsdl_s(self.Yprev, self.Y)) / sn
            epri = scipy.sqrt(self.Nc)*self.opt['AbsStopTol']/rn + \
                self.opt['RelStopTol']
            edua = scipy.sqrt(self.Nx)*self.opt['AbsStopTol']/sn + \
                self.opt['RelStopTol']

        return r, s, epri, edua 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import sqrt [as 别名]
def euclidean(h1, h2): # 9 us @array, 33 us @list \w 100 bins
    r"""
    Equal to Minowski distance with :math:`p=2`.
    
    See also
    --------
    minowski
    """
    h1, h2 = __prepare_histogram(h1, h2)
    return math.sqrt(scipy.sum(scipy.square(scipy.absolute(h1 - h2)))) 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import sqrt [as 别名]
def quadratic_forms(h1, h2):
    r"""
    Quadrativ forms metric.
    
    Notes
    -----
    UNDER DEVELOPMENT
    
    This distance measure shows very strange behaviour. The expression
    transpose(h1-h2) * A * (h1-h2) yields egative values that can not be processed by the
    square root. Some examples::
    
        h1        h2                                          transpose(h1-h2) * A * (h1-h2)
        [1, 0] to [0.0, 1.0] :                                -2.0
        [1, 0] to [0.5, 0.5] :                                 0.0
        [1, 0] to [0.6666666666666667, 0.3333333333333333] :   0.111111111111
        [1, 0] to [0.75, 0.25] :                               0.0833333333333
        [1, 0] to [0.8, 0.2] :                                 0.06
        [1, 0] to [0.8333333333333334, 0.16666666666666666] :  0.0444444444444
        [1, 0] to [0.8571428571428572, 0.14285714285714285] :  0.0340136054422
        [1, 0] to [0.875, 0.125] :                             0.0267857142857
        [1, 0] to [0.8888888888888888, 0.1111111111111111] :   0.0216049382716
        [1, 0] to [0.9, 0.1] :                                 0.0177777777778
        [1, 0] to [1, 0]:                                      0.0
    
    It is clearly undesireable to recieve negative values and even worse to get a value
    of zero for other cases than the same histograms.
    """
    h1, h2 = __prepare_histogram(h1, h2)
    A = __quadratic_forms_matrix_euclidean(h1, h2)
    return math.sqrt((h1-h2).dot(A.dot(h1-h2))) # transpose(h1-h2) * A * (h1-h2) 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import sqrt [as 别名]
def adaptive_trapezint(f, a, b, eps=1E-4):
    ddf = []
    for i in range(101):
        ddf.append(abs(diff2(f, a + i * (b - a) / 100)))
    max_ddf = max(ddf)
    h = sqrt(12 * eps) * 1 / sqrt((b - a) * max_ddf)
    n = (b - a) / h
    s = (f(a) + f(b)) / 2
    for i in range(1, int(n)):
        s += f(a + i * h)
    s *= h
    return s 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import sqrt [as 别名]
def _Efficiency(self):
        Gas = self.kwargs["entradaGas"]
        Liquido = self.kwargs["entradaLiquido"]
        rhoS = Gas.solido.rho
        muG = Gas.Gas.mu
        rhoL = Liquido.Liquido.rho

        rendimiento_fraccional = []
        if self.kwargs["modelo_rendimiento"] == 0:
            # Modelo de Johnstone (1954)
            l = sqrt(pi/8)*Gas.Gas.mu/0.4987445/sqrt(Gas.Gas.rho*Gas.P)
            for dp in Gas.solido.diametros:
                Kn = l/dp*2
                C = Cunningham(l, Kn)
                kp = C*rhoS*dp**2*self.Vg/9/Gas.Gas.mu/self.dd
                penetration = exp(-self.k*self.R*kp**0.5)
                rendimiento_fraccional.append(1-penetration)

        elif self.kwargs["modelo_rendimiento"] == 1:
            # Modelo de Calvert (1972)
            l = sqrt(pi/8)*muG/0.4987445/sqrt(Gas.Gas.rho*Gas.P)
            for dp in Gas.solido.diametros:
                Kn = l/dp*2
                C = Cunningham(l, Kn)
                kp = C*rhoS*dp**2*self.Vg/9/muG/self.dd
                b = (-0.7-kp*self.f+1.4*log((kp*self.f+0.7)/0.7)+0.49 /
                     (0.7+kp*self.f))
                penetration = exp(self.R*self.Vg*rhoL*self.dd/55/muG*b/kp)
                if penetration > 1:
                    penetration = 1
                elif penetration < 0:
                    penetration = 0
                rendimiento_fraccional.append(1-penetration)

        return rendimiento_fraccional