python类sqrt()的实例源码

def writeCgxCFile(self,filename,base_station):
        """
        write cgx type c file (without header)
        need to check which is corrected from tides
        re-write the sd as former SD/sqrt(duration)
        last column is the first gravity value of the 'o' file
        hence the first of the dictionnary, corrected for the tides

        should add a checkbox in the tree object to select the base station
        """

        offset=self.station_dic[base_station].grav[0]
        print "filename: %s"%filename
        file=open(filename,'w')

        for station_id,sta in sorted(self.station_dic.iteritems(), key=lambda x: x[1].t[1]):                  

        #for station_id,sta in self.station_dic.iteritems():
            if sta.keepitem==1:
                for i in range(len(sta.t)):
                    if sta.keepdata[i]==1:                
                        file.write("%6d %4.3f %1.3f  %2d   %2d  %2.1f   %2.1f   %2.2f   %1.3f   %3d %3.3f %02d%02d%02d %02d%02d%02d   0  0.000     %3.4f\n"%(sta.station[i],\
                            sta.grav[i]-sta.etc[i],sta.sd[i]/sqrt(sta.dur[i]),
                            sta.dur[i], sta.rej[i], sta.tiltx[i], sta.tilty[i],
                            sta.temp[i],sta.etc[i],sta.t[i].timetuple()[7],
                            sta.t[i].hour*60+sta.t[i].minute+sta.t[i].second/60,
                            sta.t[i].day,sta.t[i].month,sta.t[i].year-2000,
                            sta.t[i].hour,sta.t[i].minute,sta.t[i].second,
                              sta.grav[i]-offset))

        file.close()

def writeModifCFile(self,filename,base_station):
        """
        write modified cgx type c file (without header)
        need to check which is corrected from tides
        re-write the sd as former SD/sqrt(duration)
        last column is the first gravity value of the 'o' file
        hence the first of the dictionnary, corrected for the tides
        add an extra column with the keepdata value (0 or 1) for further 
        loading        

        should add a checkbox in the tree object to select the base station
        IMPORTANT: grav is corrected from tides (as in the classical raw data), 
        which is different than for classical c files
        and SD is SD, not SD/sqrt(dur) as for classical c files
        OR an other option could be to really write 'c' files and modify the 
        reading function (read_c_file)
        """

        offset=self.station_dic[base_station].grav[0]
        print "filename: %s"%filename
        file=open(filename,'w')

        for station_id,sta in sorted(self.station_dic.iteritems(), key=lambda x: x[1].t[1]):                  

        #for station_id,sta in self.station_dic.iteritems():
            if sta.keepitem==1:
                for i in range(len(sta.t)):               
                     file.write("%6d %4.3f %1.3f  %2d   %2d %2.1f   %2.1f   %2.2f   %1.3f   %3d %3.3f %02d%02d%02d %02d%02d%02d   0  0.000     %3.4f    %1d\n"%(sta.station[i],\
                        sta.grav[i],sta.sd[i],
                        sta.dur[i], sta.rej[i], sta.tiltx[i], sta.tilty[i],
                        sta.temp[i],sta.etc[i],sta.t[i].timetuple()[7],
                        sta.t[i].hour*60+sta.t[i].minute+sta.t[i].second/60,
                        sta.t[i].day,sta.t[i].month,sta.t[i].year-2000,
                        sta.t[i].hour,sta.t[i].minute,sta.t[i].second,
                          sta.grav[i]-offset, sta.keepdata[i]))

        file.close()

def rmse(y_test, y):  
    return sp.sqrt(sp.mean((y_test - y) ** 2))

def kde(data, N=None, MIN=None, MAX=None):

    # Parameters to set up the mesh on which to calculate
    N = 2**14 if N is None else int(2**sci.ceil(sci.log2(N)))
    if MIN is None or MAX is None:
        minimum = min(data)
        maximum = max(data)
        Range = maximum - minimum
        MIN = minimum - Range/10 if MIN is None else MIN
        MAX = maximum + Range/10 if MAX is None else MAX

    # Range of the data
    R = MAX-MIN

    # Histogram the data to get a crude first approximation of the density
    M = len(data)
    DataHist, bins = sci.histogram(data, bins=N, range=(MIN,MAX))
    DataHist = DataHist/M
    DCTData = scipy.fftpack.dct(DataHist, norm=None)

    I = [iN*iN for iN in xrange(1, N)]
    SqDCTData = (DCTData[1:]/2)**2

    # The fixed point calculation finds the bandwidth = t_star
    guess = 0.1
    try:
        t_star = scipy.optimize.brentq(fixed_point, 0, guess, 
                                       args=(M, I, SqDCTData))
    except ValueError:
        print 'Oops!'
        return None

    # Smooth the DCTransformed data using t_star
    SmDCTData = DCTData*sci.exp(-sci.arange(N)**2*sci.pi**2*t_star/2)
    # Inverse DCT to get density
    density = scipy.fftpack.idct(SmDCTData, norm=None)*N/R
    mesh = [(bins[i]+bins[i+1])/2 for i in xrange(N)]
    bandwidth = sci.sqrt(t_star)*R

    density = density/sci.trapz(density, mesh)
    return bandwidth, mesh, density

def fixed_point(t, M, I, a2):
    l=7
    I = sci.float128(I)
    M = sci.float128(M)
    a2 = sci.float128(a2)
    f = 2*sci.pi**(2*l)*sci.sum(I**l*a2*sci.exp(-I*sci.pi**2*t))
    for s in range(l, 1, -1):
        K0 = sci.prod(xrange(1, 2*s, 2))/sci.sqrt(2*sci.pi)
        const = (1 + (1/2)**(s + 1/2))/3
        time=(2*const*K0/M/f)**(2/(3+2*s))
        f=2*sci.pi**(2*s)*sci.sum(I**s*a2*sci.exp(-I*sci.pi**2*time))
    return t-(2*M*sci.sqrt(sci.pi)*f)**(-2/5)

def __init__(self, data, mleValue, fitParameters=True, mean=None, sigma=None):
        super(normalModel, self).__init__(data)




        self.MLE = mleValue

        if(None in [mean, sigma]):
            fitParameters=True
        if(fitParameters):
            mean = Mean(self.getDataSet())      
            sigma = standardDeviation(self.getDataSet())
            try:

                def normDist(x, x0, sigma):
                    output = 1.0/sqrt(2*np.pi*(sigma**2))
                    output *= exp(-0.5*((x - x0)/sigma)**2)
                    return output

                self.n, self.bins, patches = plt.hist(self.getDataSet(), self.getDatasetSize()/10, normed=1, facecolor='blue', alpha = 0.55)
                popt,pcov = curve_fit(normDist,self.bins[:-1], self.n, p0=[mean, sigma])
                ##plt.plot(bins[:-1], gaus(bins[:-1],*popt),'c-',label="Gaussian Curve with params\na=%f\nx0=%f\nsigma=%f" % (popt[0], popt[1], popt[2]), alpha=0.5)
                print "Fitted gaussian curve to data with params x0 %f, sigma %f" % (popt[0], popt[1])
                self.x0 = popt[0]
                self.sigma = popt[1]

                self.fitted = True
            except RuntimeError:
                print "Unable to fit data to normal curve, runtime error"
                raise
            except Warning:
                raise RuntimeError
        else:
            self.x0 = mean
            self.sigma = sigma

def getModelpdf(self, x):
        output = 1.0/math.sqrt(2*np.pi*(self.getSigmaValue()**2))
        output *= math.exp(-0.5*((x - self.getx0Value())/self.getSigmaValue())**2)
        return output

def getModelpdf(self, x):
        if(x>=0):
            output = 1.0/(  (x*self.getSigmaValue())*sqrt(2*np.pi) )
            output *= exp(-0.5*((log(x)- self.getx0Value())/self.getSigmaValue())**2)
        else:
            output = x      

        return scipy.where((x<0), 0.0, output)

def _calculate_whitening_transformation_matrix(self, sample_from_prior):
        samples_vec = sp.asarray([self._dict_to_to_vect(x)
                                  for x in sample_from_prior])
        # samples_vec is an array of shape nr_samples x nr_features
        means = samples_vec.mean(axis=0)
        centered = samples_vec - means
        covariance = centered.T.dot(centered)
        w, v = la.eigh(covariance)
        self._whitening_transformation_matrix = (
            v.dot(sp.diag(1. / sp.sqrt(w))).dot(v.T))

def weighted_std(points, weights):
    mean = weighted_mean(points, weights)
    std = sp.sqrt(((points - mean)**2 * weights).sum())
    return std

def test_gaussian_multiple_populations(db_path, sampler):
    sigma_x = 1
    sigma_y = .5
    y_observed = 2

    def model(args):
        return {"y": st.norm(args['x'], sigma_y).rvs()}

    models = [model]
    models = list(map(SimpleModel, models))
    nr_populations = 4
    population_size = ConstantPopulationSize(600)
    parameter_given_model_prior_distribution = [Distribution(x=st.norm(0, sigma_x))]
    abc = ABCSMC(models, parameter_given_model_prior_distribution,
                 MinMaxDistanceFunction(measures_to_use=["y"]),
                 population_size,
                 eps=MedianEpsilon(.2),
                 sampler=sampler)

    abc.new(db_path, {"y": y_observed})

    minimum_epsilon = -1

    abc.do_not_stop_when_only_single_model_alive()
    history = abc.run(minimum_epsilon, max_nr_populations=nr_populations)
    posterior_x, posterior_weight = history.get_distribution(0, None)
    posterior_x = posterior_x["x"].as_matrix()
    sort_indices = sp.argsort(posterior_x)
    f_empirical = sp.interpolate.interp1d(sp.hstack((-200, posterior_x[sort_indices], 200)),
                                          sp.hstack((0, sp.cumsum(posterior_weight[sort_indices]), 1)))

    sigma_x_given_y = 1 / sp.sqrt(1 / sigma_x**2 + 1 / sigma_y**2)
    mu_x_given_y = sigma_x_given_y**2 * y_observed / sigma_y**2
    expected_posterior_x = st.norm(mu_x_given_y, sigma_x_given_y)
    x = sp.linspace(-8, 8)
    max_distribution_difference = sp.absolute(f_empirical(x) - expected_posterior_x.cdf(x)).max()
    assert max_distribution_difference < 0.052
    assert history.max_t == nr_populations-1
    mean_emp, std_emp = mean_and_std(posterior_x, posterior_weight)
    assert abs(mean_emp - mu_x_given_y) < .07
    assert abs(std_emp - sigma_x_given_y) < .12

def test_two_competing_gaussians_single_population(db_path, sampler, transition):
    sigma_x = .5
    sigma_y = .5
    y_observed = 1

    def model(args):
        return {"y": st.norm(args['x'], sigma_y).rvs()}

    models = [model, model]
    models = list(map(SimpleModel, models))
    population_size = ConstantPopulationSize(500)
    mu_x_1, mu_x_2 = 0, 1
    parameter_given_model_prior_distribution = [
        Distribution(x=st.norm(mu_x_1, sigma_x)),
        Distribution(x=st.norm(mu_x_2, sigma_x))]
    abc = ABCSMC(models, parameter_given_model_prior_distribution,
                 MinMaxDistanceFunction(measures_to_use=["y"]),
                 population_size,
                 transitions=[transition(), transition()],
                 eps=MedianEpsilon(.02),
                 sampler=sampler)
    abc.new(db_path, {"y": y_observed})

    minimum_epsilon = -1
    nr_populations = 1
    abc.do_not_stop_when_only_single_model_alive()
    history = abc.run(minimum_epsilon, max_nr_populations=1)
    mp = history.get_model_probabilities(history.max_t)

    def p_y_given_model(mu_x_model):
        return st.norm(mu_x_model, sp.sqrt(sigma_y**2+sigma_x**2)).pdf(y_observed)

    p1_expected_unnormalized = p_y_given_model(mu_x_1)
    p2_expected_unnormalized = p_y_given_model(mu_x_2)
    p1_expected = p1_expected_unnormalized / (p1_expected_unnormalized + p2_expected_unnormalized)
    p2_expected = p2_expected_unnormalized / (p1_expected_unnormalized + p2_expected_unnormalized)
    assert history.max_t == nr_populations - 1
    assert abs(mp.p[0] - p1_expected) + abs(mp.p[1] - p2_expected) < .07

def rmse(y_test, y):
    return sp.sqrt(sp.mean((y_test - y) ** 2))

def distance(p1, p2):
    return sqrt((p1[0]-p2[0])**2+(p1[1]-p2[1])**2+(p1[2]-p2[2])**2)

def bsCall(S,X,T,r,sigma):
    from scipy import log,exp,sqrt,stats
    d1=(log(S/X)+(r+sigma*sigma/2.)*T)/(sigma*sqrt(T))
    d2 = d1-sigma*sqrt(T)
    return S*stats.norm.cdf(d1)-X*exp(-r*T)*stats.norm.cdf(d2)

def sharpe(R,w):
    var = portfolio_var(R,w)
    mean_return=sp.mean(R,axis=0)
    ret = sp.array(mean_return)
    return (sp.dot(w,ret) - rf)/sp.sqrt(var)
# function 4: for given n-1 weights, return a negative sharpe ratio

def BIS_f(R,s,n):
    R=R0
    for i in sp.arange(0,n):
        deltaR=z[i]*s/sp.sqrt(2.)
        logR=sp.log(R)
        R=sp.exp(logR+deltaR)
        output.append(round(R,5))
    return output 
#

def bsCall(S,X,T,r,sigma):
    d1=(log(S/X)+(r+sigma*sigma/2.)*T)/(sigma*sqrt(T)) 
    d2 = d1-sigma*sqrt(T)
    return S*stats.norm.cdf(d1)-X*exp(-r*T)*stats.norm.cdf(d2)
#

def bs_call(S,X,T,r,sigma):
    d1=(log(S/X)+(r+sigma*sigma/2.)*T)/(sigma*sqrt(T)) 
    d2 = d1-sigma*sqrt(T)
    return S*stats.norm.cdf(d1)-X*exp(-r*T)*stats.norm.cdf(d2)

def implied_vol_call(S,X,T,r,c):
    for i in range(200):
        sigma=0.005*(i+1)
        d1=(log(S/X)+(r+sigma*sigma/2.)*T)/(sigma*sqrt(T))
        d2 = d1-sigma*sqrt(T)
        diff=c-(S*stats.norm.cdf(d1)-X*exp(-r*T)*stats.norm.cdf(d2))
        if abs(diff)<=0.01:
            return i,sigma, diff