Python pylab.std函数代码示例

def scatter_stats(db, s1, s2, f1=None, f2=None, **kwargs):
    if f1 == None:
        f1 = lambda x: x  # constant function

    if f2 == None:
        f2 = f1

    x = []
    xerr = []

    y = []
    yerr = []

    for k in db:
        x_k = [f1(x_ki) for x_ki in db[k].__getattribute__(s1).gettrace()]
        y_k = [f2(y_ki) for y_ki in db[k].__getattribute__(s2).gettrace()]

        x.append(pl.mean(x_k))
        xerr.append(pl.std(x_k))

        y.append(pl.mean(y_k))
        yerr.append(pl.std(y_k))

        pl.text(x[-1], y[-1], " %s" % k, fontsize=8, alpha=0.4, zorder=-1)

    default_args = {"fmt": "o", "ms": 10}
    default_args.update(kwargs)
    pl.errorbar(x, y, xerr=xerr, yerr=yerr, **default_args)
    pl.xlabel(s1)
    pl.ylabel(s2)

def data_to_ch(data):
    ch = {}
    for ch_ind in range(1, 97):
        ch[ch_ind] = {}
        ch[ch_ind]["bl"] = data[ch_ind]["blanks"]
        ch[ch_ind]["bl_mu"] = pl.mean(ch[ch_ind]["bl"])
        ch[ch_ind]["bl_sem"] = pl.std(ch[ch_ind]["bl"]) / pl.sqrt(len(ch[ch_ind]["bl"]))
        for ind in sorted(data[ch_ind].keys()):
            if ind != "blanks":
                k = ind[0]
                if k not in ch[ch_ind]:
                    ch[ch_ind][k] = {}
                    ch[ch_ind][k]["fr"] = []
                    ch[ch_ind][k]["fr_mu"] = []
                    ch[ch_ind][k]["fr_sem"] = []
                    ch[ch_ind][k]["pos_y"] = []
                    ch[ch_ind][k]["dprime"] = []
                ch[ch_ind][k]["fr"].append(data[ch_ind][ind]["on"])
                ch[ch_ind][k]["fr_mu"].append(pl.mean(data[ch_ind][ind]["on"]))
                ch[ch_ind][k]["fr_sem"].append(pl.std(data[ch_ind][ind]["on"]) / pl.sqrt(len(data[1][ind]["on"])))
                ch[ch_ind][k]["pos_y"].append(ind[2])
                # print ch[ch_ind][k]['pos_y']
                # print pl.std(data[ch_ind][ind]['on'])
                ch[ch_ind][k]["dprime"].append(
                    (pl.mean(data[ch_ind][ind]["on"]) - ch[ch_ind]["bl_mu"])
                    / ((pl.std(ch[ch_ind]["bl"]) + pl.std(data[ch_ind][ind]["on"])) / 2)
                )
                # print ch[ch_ind]['OSImage_5']['pos_y']
    return ch

def compare_models(db, stoch="itn coverage", stat_func=None, plot_type="", **kwargs):
    if stat_func == None:
        stat_func = lambda x: x

    X = {}
    for k in sorted(db.keys()):
        c = k.split("_")[2]
        X[c] = []

    for k in sorted(db.keys()):
        c = k.split("_")[2]
        X[c].append([stat_func(x_ki) for x_ki in db[k].__getattribute__(stoch).gettrace()])

    x = pl.array([pl.mean(xc[0]) for xc in X.values()])
    xerr = pl.array([pl.std(xc[0]) for xc in X.values()])
    y = pl.array([pl.mean(xc[1]) for xc in X.values()])
    yerr = pl.array([pl.std(xc[1]) for xc in X.values()])

    if plot_type == "scatter":
        default_args = {"fmt": "o", "ms": 10}
        default_args.update(kwargs)
        for c in X.keys():
            pl.text(pl.mean(X[c][0]), pl.mean(X[c][1]), " %s" % c, fontsize=8, alpha=0.4, zorder=-1)
        pl.errorbar(x, y, xerr=xerr, yerr=yerr, **default_args)
        pl.xlabel("First Model")
        pl.ylabel("Second Model")
        pl.plot([0, 1], [0, 1], alpha=0.5, linestyle="--", color="k", linewidth=2)

    elif plot_type == "rel_diff":
        d1 = sorted(100 * (x - y) / x)
        d2 = sorted(100 * (xerr - yerr) / xerr)
        pl.subplot(2, 1, 1)
        pl.title("Percent Model 2 deviates from Model 1")

        pl.plot(d1, "o")
        pl.xlabel("Countries sorted by deviation in mean")
        pl.ylabel("deviation in mean (%)")

        pl.subplot(2, 1, 2)
        pl.plot(d2, "o")
        pl.xlabel("Countries sorted by deviation in std err")
        pl.ylabel("deviation in std err (%)")
    elif plot_type == "abs_diff":
        d1 = sorted(x - y)
        d2 = sorted(xerr - yerr)
        pl.subplot(2, 1, 1)
        pl.title("Percent Model 2 deviates from Model 1")

        pl.plot(d1, "o")
        pl.xlabel("Countries sorted by deviation in mean")
        pl.ylabel("deviation in mean")

        pl.subplot(2, 1, 2)
        pl.plot(d2, "o")
        pl.xlabel("Countries sorted by deviation in std err")
        pl.ylabel("deviation in std err")
    else:
        assert 0, "plot_type must be abs_diff, rel_diff, or scatter"

    return pl.array([x, y, xerr, yerr])

def plot2():
    import pylab as pl
    hs, ds = [], []
    for event, time in load():
        if event == main_start:
            start_time = time
        elif event == main_end:
            d0, h0 = days_hours(start_time)
            d1, h1 = days_hours(time)
            hs.append((h0, h1))
            ds.append((d0, d1))
            pl.plot([d0, d1], [h0, h1], 'b')
    ihs, fhs = zip(*hs)
    ids, fds = zip(*ds)
    pl.plot(ids, ihs, 'g')
    pl.plot([ids[0], ids[-1]], [pl.mean(ihs)] * 2, 'g--')
    pl.plot(fds, fhs, 'r')
    pl.plot([fds[0], fds[-1]], [pl.mean(fhs)] * 2, 'r--')
    f, i = pl.mean(fhs), pl.mean(ihs)
    pl.plot([fds[0], fds[-1]], [(f + i) / 2] * 2, 'b--')
    print i, f, f - i, (f + i) / 2
    std_i, std_f = pl.std(ihs), pl.std(fhs)
    print std_i, std_f
    pl.xlim(ids[0], fds[-1])
    pl.ylim(4, 28)
    pl.grid(True)
    pl.xlabel('Time [day]')
    pl.ylabel('Day interval [hours]')
    pl.show()

 def stderr(X,Y=None):
     if len(X) <= 1: return 0.0
     stderr_x = pow(pylab.std(X),2)/len(X)
     if Y:
         if len(Y) <= 1: return 0.0
         stderr_y = pow(pylab.std(Y),2)/len(Y)
     else: stderr_y = 0
     return math.sqrt(stderr_x + stderr_y)

def Vplot(Ws):
	"""Calculate the potential function and plot it"""

	N_bstrp = input("Please enter the number of bootstraps: ")
	N_bin = input("Please enter the bin size: ")
	style = raw_input("Please enter a linestyle: ")

	Ws = bin(Ws,N_bin)
	aVs = pl.zeros((N_bstrp,) + pl.shape(Ws)[1:])
	bs = pl.zeros((N_bstrp,3))
        
	for i in xrange(N_bstrp):
		W = pl.mean(bootstrap(Ws),axis=0)
		aVs[i] = calcaV(W,method="fit")
		bs[i] = potfit(aVs[i,:,0])

			
	r = pl.arange(1,7)
	aV = pl.mean(aVs,axis=0)
	aVerr = pl.std(aVs,axis=0)
	b = pl.mean(bs,axis=0)

	a_s = 0.5 / pl.sqrt((1.65 + bs[:,1]) / bs[:,0])
	sigmas = bs[:,0] / a_s**2
	Bs = bs[:,1]
	As = bs[:,2] / a_s

	a = pl.mean(a_s)
	aerr = pl.std(a_s)
	sigma = pl.mean(sigmas)
	sigmaerr = pl.std(sigmas)
	B = pl.mean(Bs)
	Berr = pl.std(Bs)
	A = pl.mean(As)
	Aerr = pl.std(As)

	print("Fit parameters:")
	print("sigma = %f +/- %f fm^-2 = %f +/- %f MeV^2"
		% (sigma, sigmaerr, sigma * 197**2, sigmaerr * 197**2))
	print("B = %f +/- %f" % (B, Berr))
	print("A = %f +/- %f fm^-1 = %f +/- %f MeV"
		% (A, Aerr, A*197, Aerr*197))
	print("Lattice spacing, a = %f +/- %f fm = %f +/- %f MeV^-1"
		% (a, aerr, a/197, aerr/197))
	
	r_fit = pl.arange(0.25,r[-1]+1,0.1)
	aV_fit = V(b,r_fit)
	
	handles = []
	handles.append(pl.errorbar(r,aV[:,0],yerr=aVerr[:,0],fmt='o'+style[0]))
	handles.append(pl.plot(r_fit,aV_fit,style))
	pl.ylim([0,pl.nanmax(aV)+0.25])
	pl.xlim([0,pl.nanmax(r_fit)+0.25])
	pl.xlabel("$r / a$")
	pl.ylabel("$aV(r)$")

	return aV,handles

 def update(self,t,val):
     oldavg = self.avg
     avrg = Ema.update(self,t,val)
     
     self.__samples.append((t,self.lastvalue))
     self.__samples_nonone.append(self.lastvalue)
     
     if oldavg == None:
         self.first_t = t 
         return (None,None,None)
     newavg = avrg
     
     # check limits of timeframe 
     while t - self.__samples[0][0] > self.timeframe and len(self.__samples) > 2:
         _, pv = self.__samples.pop(0)
         del self.__samples_nonone[0]
         #self.variance += (val-pv)*(val-newavg+pv-oldavg)/(self.timeframe) # this seems to use constant number of samples
         #std = math.sqrt(self.variance)
         
     std = pylab.std(self.__samples_nonone)*self.k # this takes long
     
     if avrg != None:
         return avrg,avrg+std,avrg-std
     else:
         return (None,None,None)

def flow_rate_hist(sheets):
    ant_rates = []
    weights = []
    for sheet in sheets:
        ants, seconds, weight = flow_rate(sheet)
        ant_rate = seconds / ants
        #ant_rate = ants / seconds
        ant_rates.append(ant_rate)
        weights.append(float(weight))
        #weights.append(seconds)

    weights = pylab.array(weights)
    weights /= sum(weights)

    #print "ants per second"
    print "seconds per ant"
    mu = pylab.mean(ant_rates)
    print "mean", pylab.mean(ant_rates)
    wmean = pylab.average(ant_rates, weights=weights)
    print "weighted mean", wmean
    print "median", pylab.median(ant_rates)
    print "std", pylab.std(ant_rates, ddof=1)
    ant_rates = pylab.array(ant_rates)
    werror = (ant_rates - mu) * weights
    print "weighted std", ((sum(werror ** 2))) ** 0.5
    print "weighted std 2", (pylab.average((ant_rates - mu)**2, weights=weights)) ** 0.5
    pylab.figure()
    pylab.hist(ant_rates)
    pylab.savefig('ant_flow_rates.pdf', format='pdf')
    pylab.close()

def _CalcMutualNearestNeighbors(hull_points, all_points):
    all_points_list = list(all_points)
    ds = distance.pdist(list(all_points))
    std_d = p.std(ds)
    
    square_ds = distance.squareform(ds)
    nearest_neighbors = {}
    
    for i, point in enumerate(all_points_list):
        if point not in hull_points:
            continue
        
        my_ds = [(d, j) for j, d in enumerate(square_ds[i])
                 if j != i]
        my_ds.sort()
        nearest_neighbors[point] = set([j for d,j in my_ds[:3]])
    
    no_mutual = set()
    for i, point in enumerate(all_points_list):
        if point not in hull_points:
            continue
        
        no_nbrs = True
        for neighbor_index in nearest_neighbors.get(point, []):
            neighbor = all_points_list[neighbor_index]
            neighbor_set = nearest_neighbors.get(neighbor, [])
            if i in neighbor_set:
                no_nbrs = False
        
        if no_nbrs:
            no_mutual.add(point)
                
    return no_mutual

 def generate_normalized_test_data(self,
                                   channels, 
                                   time_points, 
                                   function, 
                                   sampling_frequency, 
                                   initial_phase=0.0):
     """
     A method which generates a normalized (mu = 0, sigma =1) signal for testing, with
     the specified number of "channels" which are all generated using the given function
     """
    
     #Generate an empty ndarray
     data = numpy.zeros((time_points, channels))
     
     #Compute the values for all channels
     for channel_index in range(channels):
         for time_index in range(time_points):
             data[time_index, channel_index] = function(2.0 * numpy.pi * (channel_index + 1) * (time_index / sampling_frequency + initial_phase))
         current_channel = data[:, channel_index]
         current_channel = (current_channel - pylab.mean(current_channel))/pylab.std(current_channel)
         data[:, channel_index] = current_channel
         
     #Generate a time series build out of the data
     test_data = TimeSeries(input_array = data, 
                            channel_names = [("test_channel_%s" % i) for i in range(channels)],
                            sampling_frequency = sampling_frequency,
                            start_time = initial_phase,
                            end_time = float(time_points) / sampling_frequency + initial_phase)
     
     return test_data

def pc_pm_std(data, ndim):
    """
    This is a helper function.
    It returns the value of +1 * std(x), where x is the ndim-th principal
    component of the data

    Parameters:
    -----------
    data: `array` (*n*-by-*d*)
        the data on which the principal component analysis is performed.
    ndim: `integer`
        the number of the principal axis on which the analysis is performed.
        **NOTE** this is zero-based, i.e. to compute the first principal
        component, ndim=0

    Returns:
    --------
    std_pc: `array` (1-by-*d*)
        the vector that points in the direction of the *ndim*th principal
        axis, and has the length of the standard deviation of the scores
        along this axis.

    """

    u,s,v = svd(data.T, full_matrices = False)
    direction = u[:, ndim : ndim + 1]
    scale = std(dot(direction.T, data.T))
    return scale * direction.T

    def zeroPaddData(self,desiredLength,paddmode='zero',where='end'):    
        #zero padds the time domain data, it is possible to padd at the beginning,
        #or at the end, and further gaussian or real zero padding is possible        
        #might not work for gaussian mode!

        desiredLength=int(desiredLength)
        #escape the function        
        if desiredLength<0:
            return 0

        #calculate the paddvectors        
        if paddmode=='gaussian':
            paddvec=py.normal(0,py.std(self.getPreceedingNoise())*0.05,desiredLength)
        else:
            paddvec=py.ones((desiredLength,self.tdData.shape[1]-1))
            paddvec*=py.mean(self.tdData[-20:,1:])
            
        timevec=self.getTimes()
        if where=='end':
            #timeaxis:
            newtimes=py.linspace(timevec[-1],timevec[-1]+desiredLength*self.dt,desiredLength)
            paddvec=py.column_stack((newtimes,paddvec))
            longvec=py.row_stack((self.tdData,paddvec))
        else:
            newtimes=py.linspace(timevec[0]-(desiredLength+1)*self.dt,timevec[0],desiredLength)
            paddvec=py.column_stack((newtimes,paddvec))
            longvec=py.row_stack((paddvec,self.tdData))
            
        self.setTDData(longvec)

 def post_lecture(self):
     STD = std(self.Y, 1)
     MM = mean(self.Y, 1)
     TT, self.NN = self.Y.shape
     if self.centred:
         for t in xrange(0, TT):
             self.Y[t, :] = (self.Y[t, :] - MM[t]) / STD[t]

 def getelnNoise(self,tdData):
     #returns the uncertainty due to electronic noise
     
     #signal preceeding the pulse (X and Y channel)  
     precNoise=self.getPreceedingNoise(tdData)
     #is this normalization really correct?!
     elnNoise = py.std(precNoise, ddof = 1,axis=0)/py.sqrt(precNoise.shape[0])
     return elnNoise

def readDatDirectory(key, directory):
    global stats
    #Don't read data in if it's already read
    if not key in DATA["mean"]:
        data = defaultdict(array)

        #Process the dat files
        for datfile in glob.glob(directory + "/*.dat"):
            fileHandle = open(datfile, 'rb')
            keys, dataDict = csvExtractAllCols(fileHandle)
            stats = union(stats, keys)
            for aKey in keys:
                if not aKey in data:
                    data[aKey] = reshape(array(dataDict[aKey]),
                                         (1, len(dataDict[aKey])))
                else:
                    data[aKey] = append(data[aKey],
                                        reshape(array(dataDict[aKey]),
                                                (1, len(dataDict[aKey]))),
                                        axis=0)

        #Process the div files'
        for datfile in glob.glob(directory + "/*.div"):
            fileHandle = open(datfile, 'rb')
            keys, dataDict = csvExtractAllCols(fileHandle)
            stats = union(stats, keys)
            for aKey in keys:
                if not aKey in data:
                    data[aKey] = reshape(array(dataDict[aKey]),
                                         (1, len(dataDict[aKey])))
                else:
                    data[aKey] = append(data[aKey],
                                        reshape(array(dataDict[aKey]),
                                                (1, len(dataDict[aKey]))),
                                        axis=0)

        #Iterate through the stats and calculate mean/standard deviation
        for aKey in stats:
            if aKey in data:
                DATA["mean"][key][aKey] = mean(data[aKey], axis=0)
                DATA["median"][key][aKey] = median(data[aKey], axis=0)
                DATA["std"][key][aKey] = std(data[aKey], axis=0)
                DATA["ste"][key][aKey] = std(data[aKey], axis=0)/ sqrt(len(data[aKey]))
                DATA["min"][key][aKey] = mean(data[aKey], axis=0)-amin(data[aKey], axis=0)
                DATA["max"][key][aKey] = amax(data[aKey], axis=0)-mean(data[aKey], axis=0)
                DATA["actual"][key][aKey] = data[aKey]