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
python类sqrt()的实例源码
def lsInversion(self):
"""
LS Inversion from Hwang et al (2002)
"""
At=np.transpose(self.A)
St=np.transpose(self.S)
N=At.dot(self.P).dot(self.A)
#solution:
self.X=np.linalg.inv(N+self.S.dot(St)).dot(At).dot(self.P).dot(self.Obs)
self.r=self.A.dot(self.X)-self.Obs
rt=np.transpose(self.r)
self.VtPV=rt.dot(self.P).dot(self.r)
var_post_norm=self.VtPV/self.dof
self.SDaposteriori=np.sqrt(var_post_norm)
cov_post=np.linalg.inv(N)*var_post_norm
self.var=np.diagonal(cov_post)
def lsStatistics(self,alpha):
"""
a priori variance of unit weight = 1
"""
self.chi2=self.VtPV
t=np.sqrt(2*np.log(1/alpha))
chi_1_alpha=t-(2.515517+0.802853*t+0.010328*t**2)/(1+1.432788*t+0.189269*t**2+0.001308*t**3)
dof=float(self.dof)
self.chi2c=dof*(chi_1_alpha*sqrt(2/(9*dof))+1-2/(9*dof))**3
print "SD a posteriori: %6.2f"%float(self.SDaposteriori)
print "chi square value: %6.2f"%float(self.chi2)
print "critical chi square value: %6.2f"%float(self.chi2c)
if self.chi2<self.chi2c:
print "Chi-test accepted"
else:
print "Chi-test rejected"
def loss(r_0, a, alpha, x_focus, x, y_focus, y, overlap):
loss = np.zeros(np.size(x))
loss_squared = np.zeros(np.size(x))
dx = np.zeros(len(x))
dy = np.zeros(len(y))
for i in range(0, np.size(x)):
dx = x_focus-x[i]
dy = abs(y_focus-y[i])
R = r_0[i]+dx*alpha
if dx > 0:
loss[i] = overlap[i]*2.*a*(r_0[i]/(r_0[i]+alpha*(dx)))**2*0.5*(np.cos(-dy*pi/(R+4*r_0[i])))
loss_squared[i] = loss[i]**2
else:
loss[i] = 0
loss_squared[i] = 0
total_loss = sp.sqrt(np.sum(loss_squared))
return total_loss
poisson_disk_sampling.py 文件源码
项目:nodebox1-generative-tools
作者: x-raizor
项目源码
文件源码
阅读 46
收藏 0
点赞 0
评论 0
def __init__( self, w, h, r, n ):
# w and h are the width and height of the field
self.w = w
self.h = h
# n is the number of test points
self.n = n
self.r2 = r**2.0
self.A = 3.0*self.r2
# cs is the cell size
self.cs = r / scipy.sqrt(2)
# gw and gh are the number of grid cells
self.gw = int( scipy.ceil( self.w/self.cs ) )
self.gh = int( scipy.ceil( self.h/self.cs ) )
# create a grid and a queue
self.grid = [ None ] * self.gw * self.gh
self.queue = list()
# set the queue size and sample size to zero
self.qs, self.ss = 0, 0
poisson_disk_sampling.py 文件源码
项目:nodebox1-generative-tools
作者: x-raizor
项目源码
文件源码
阅读 28
收藏 0
点赞 0
评论 0
def rvs( self ):
if self.ss == 0:
x = random.random() * self.w
y = random.random() * self.h
self.set_point( x, y )
while self.qs:
x_idx = int( random.random() * self.qs )
s = self.queue[ x_idx ]
for y_idx in range( self.n ):
a = 2 * scipy.pi * random.random()
b = scipy.sqrt( self.A * random.random() + self.r2 )
x = s[0] + b*scipy.cos( a )
y = s[1] + b*scipy.sin( a )
if( x >= 0 )and( x < self.w ):
if( y >= 0 )and( y < self.h ):
if( self.distance( x, y ) ):
self.set_point( x, y )
del self.queue[x_idx]
self.qs -= 1
sample = list( filter( None, self.grid ) )
sample = scipy.asfarray( sample )
return sample
def fit(self, X, w):
if len(X) == 0:
raise NotEnoughParticles("Fitting not possible.")
self.X_arr = X.as_matrix()
ctree = cKDTree(X)
_, indices = ctree.query(X, k=min(self.k + 1, X.shape[0]))
covs, inv_covs, dets = list(zip(*[self._cov_and_inv(n, indices)
for n in range(X.shape[0])]))
self.covs = sp.array(covs)
self.inv_covs = sp.array(inv_covs)
self.determinants = sp.array(dets)
self.normalization = sp.sqrt(
(2 * sp.pi) ** self.X_arr.shape[1] * self.determinants)
if not sp.isreal(self.normalization).all():
raise Exception("Normalization not real")
self.normalization = sp.real(self.normalization)
c13_08_KMF_function.py 文件源码
项目:Python-for-Finance-Second-Edition
作者: PacktPublishing
项目源码
文件源码
阅读 33
收藏 0
点赞 0
评论 0
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
#
c10_33_implied_vol_EuropeanPut_min.py 文件源码
项目:Python-for-Finance-Second-Edition
作者: PacktPublishing
项目源码
文件源码
阅读 39
收藏 0
点赞 0
评论 0
def implied_vol_put_min(S,X,T,r,p):
implied_vol=1.0
min_value=100.0
for i in xrange(1,10000):
sigma=0.0001*(i+1)
d1=(log(S/X)+(r+sigma*sigma/2.)*T)/(sigma*sqrt(T))
d2 = d1-sigma*sqrt(T)
put=X*exp(-r*T)*stats.norm.cdf(-d2)-S*stats.norm.cdf(-d1)
abs_diff=abs(put-p)
if abs_diff<min_value:
min_value=abs_diff
implied_vol=sigma
k=i
put_out=put
print 'k, implied_vol, put, abs_diff'
return k,implied_vol, put_out,min_value
c12_19_up_and_out_call.py 文件源码
项目:Python-for-Finance-Second-Edition
作者: PacktPublishing
项目源码
文件源码
阅读 50
收藏 0
点赞 0
评论 0
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
#
c14_17_up_call.py 文件源码
项目:Python-for-Finance-Second-Edition
作者: PacktPublishing
项目源码
文件源码
阅读 34
收藏 0
点赞 0
评论 0
def upCall(s,x,T,r,sigma,nSimulation,barrier):
import scipy as sp
import p4f
n_steps=100
dt=T/n_steps
inTotal=0
outTotal=0
for j in range(0, nSimulation):
sT=s
inStatus=False
outStatus=True
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:
outStatus=False
inStatus=True
if outStatus==True:
outTotal+=p4f.bs_call(s,x,T,r,sigma)
else:
inTotal+=p4f.bs_call(s,x,T,r,sigma)
return outTotal/nSimulation, inTotal/nSimulation
#
c14_20_lookback_min_price_as_strike.py 文件源码
项目:Python-for-Finance-Second-Edition
作者: PacktPublishing
项目源码
文件源码
阅读 42
收藏 0
点赞 0
评论 0
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
#
c14_16_up_and_out_call.py 文件源码
项目:Python-for-Finance-Second-Edition
作者: PacktPublishing
项目源码
文件源码
阅读 28
收藏 0
点赞 0
评论 0
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
#
c14_26_up_and_out_call2.py 文件源码
项目:Python-for-Finance-Second-Edition
作者: PacktPublishing
项目源码
文件源码
阅读 51
收藏 0
点赞 0
评论 0
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
#
def _BlackScholesCall(S, K, T, sigma, r, q):
d1 = (log(S/K) + (r - q + (sigma**2)/2)*T)/(sigma*sqrt(T))
d2 = d1 - sigma*sqrt(T)
return S*exp(-q*T)*norm.cdf(d1) - K*exp(-r*T)*norm.cdf(d2)
def _BlackScholesPut(S, K, T, sigma, r, q):
d1 = (log(S/K) + (r - q + (sigma**2)/2)*T)/(sigma*sqrt(T))
d2 = d1 - sigma*sqrt(T)
return K*exp(-r*T)*norm.cdf(-d2) - S*exp(-q*T)*norm.cdf(-d1)
def _fprime(self, sigma):
logSoverK = log(self.S/self.K)
n12 = ((self.r + sigma**2/2)*self.T)
numerd1 = logSoverK + n12
d1 = numerd1/(sigma*sqrt(self.T))
return self.S*sqrt(self.T)*norm.pdf(d1)*exp(-self.r*self.T)
def computeFaceRecord(self,img,rect=None,eyes=None):
'''Given an image and face detection box, compute a face identification record'''
assert self.trained
img = self.cropFace(img,eyes)
vec = self.computeVector(img)
fir = self.pca.project(vec,whiten=True)
if self.measure == PCA_COS:
scale = scipy.sqrt((fir*fir).sum())
fir = (1.0/scale)*fir
return fir
def similarity(self,fir1,fir2):
'''Compute the similarity of two faces'''
assert self.trained == True
if self.measure == PCA_L1:
return (scipy.abs(fir1-fir2)).sum()
if self.measure == PCA_L2:
return scipy.sqrt(((fir1-fir2)*(fir1-fir2)).sum())
if self.measure == PCA_COS:
return (fir1*fir2).sum()
raise NotImplementedError("Unknown distance measure: %d"%self.measure)
def norm(self):
x = centered1d(self.x)
y = centered1d(self.y)
return scipy.sqrt(trapz2(self.intensity(), x=x, y=y))
def update_rho(self, k, r, s):
"""Automatic rho adjustment."""
if self.opt['AutoRho', 'Enabled']:
tau = self.rho_tau
mu = self.rho_mu
xi = self.rho_xi
if k != 0 and scipy.mod(k+1, self.opt['AutoRho', 'Period']) == 0:
if self.opt['AutoRho', 'AutoScaling']:
if s == 0.0 or r == 0.0:
rhomlt = tau
else:
rhomlt = scipy.sqrt(r/(s*xi) if r > s*xi else (s*xi)/r)
if rhomlt > tau:
rhomlt = tau
else:
rhomlt = tau
rsf = 1.0
if r > xi*mu*s:
rsf = rhomlt
elif s > (mu/xi)*r:
rsf = 1.0/rhomlt
self.rho = self.dtype.type(rsf*self.rho)
self.U /= rsf
if rsf != 1.0:
self.rhochange()
def rsdl_rn(self, AX, Y):
"""Compute primal residual normalisation term.
Overriding this method is required if methods :meth:`cnst_A`,
:meth:`cnst_AT`, :meth:`cnst_B`, and :meth:`cnst_c` are not
overridden.
"""
if not hasattr(self, '_cnst_nrm_c'):
self._cnst_nrm_c = np.sqrt(linalg.norm(self.cnst_c0())**2 +
linalg.norm(self.cnst_c1())**2)
return max((linalg.norm(AX), linalg.norm(Y), self._cnst_nrm_c))
def rsdl_s(self, Yprev, Y):
"""Compute dual residual vector."""
# Since s = rho A^T B (y^(k+1) - y^(k)) and B = -(I I I ...)^T,
# the correct calculation here would involve replicating (Yprev - Y)
# on the axis on which the blocks of X are stacked. Since this would
# require allocating additional memory, and since it is only the norm
# of s that is required, instead of replicating the vector it is
# scaled to have the same l2 norm as the replicated vector
return np.sqrt(self.Nb)*self.rho*(Yprev - Y)
def rsdl_rn(self, AX, Y):
"""Compute primal residual normalisation term."""
# The primal residual normalisation term is
# max( ||A x^(k)||_2, ||B y^(k)||_2 ) and B = -(I I I ...)^T.
# The scaling by sqrt(Nb) of the l2 norm of Y accounts for the
# block replication introduced by multiplication by B
return max((linalg.norm(AX), np.sqrt(self.Nb)*linalg.norm(Y)))
def c1(R1):
return c2 * sc.sqrt(R2 / R1)
# Equilibrium magnetic field strength within the slab:
def m0(W):
return sc.sqrt((c0**2 - W**2)*(vA**2 - W**2) /
((c0**2 + vA**2)*(cT**2 - W**2)))
def m00(W):
return sc.sqrt(1 - W**2/c0**2)
def m1(W, R1):
return sc.sqrt(1 - W**2*(c2*sc.sqrt(R2/R1))**(-2))
def m2(W):
return sc.sqrt(1 - W**2/c2**2)
linear_regression.py 文件源码
项目:Stock-SentimentAnalysis
作者: JoshuaMichaelKing
项目源码
文件源码
阅读 45
收藏 0
点赞 0
评论 0
def rmse(ymodel, yreal):
'''
root-mean-square-error
????????????(??????????????????)
'''
return sp.sqrt(sp.mean((ymodel - yreal) ** 2))