Python pylab.mean函数代码示例

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 broadgauss(x, y, sigma):
    '''Gaussian function for broadening
    '''

    bla = True
    plot = False
    c = 299792458.

    if bla:
        print " sigma = ", round(sigma, 4), " km/s"

    sigma = sigma * 1.0e3/c * pl.mean(x)   # sigma in Å

    if bla:
        print " sigma = ", round(sigma, 3), "  Å "

    xk = x - pl.mean(x)

    g = make_gauss(1, 0, sigma)
    yk = [g(i) for i in xk]

    if bla:
        print " Integral of the gaussian function: ", pl.trapz(yk, xk).__format__('5.3')

    if plot:
        pl.figure(2)
        pl.plot(xk, yk, '+-')
        pl.show()
    #if bla: print" size y:", y.size
    y = pl.convolve(y, yk, mode='same')
    #if bla: print" size y:", y.size

    return y/max(y)

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 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 calZsocre(self,core,surface,sampleSize):
     coreMean=mean(core)
     s=[]
     for i in range(sampleSize):
         s.append(mean(sample(surface,len(core))))
     sig= sqrt(var(s))
     return (coreMean-mean(s))/sig

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 latent_simplex(X):
    """ TODO: describe this function"""
    N, T, J = X.shape

    alpha = []
    for t in range(T):
        alpha_t = []
        for j in range(J):
            mu_alpha_tj = pl.mean(X[:,t,j]) / pl.mean(X[:,t,:], 0).sum()
            alpha_t.append(mc.Normal('alpha_%d_%d'%(t,j), mu=0., tau=1., value=pl.log(mu_alpha_tj)))
        alpha.append(alpha_t)

    @mc.deterministic
    def pi(alpha=alpha):
        pi = pl.zeros((T, J))
        for t in range(T):
            pi[t] = pl.reshape(pl.exp(alpha[t]), J) / pl.sum(pl.exp(alpha[t]))
        return pi

    @mc.observed
    def X_obs(pi=pi, value=X.mean(0), sigma=X.std(0), pow=2):
        """ TODO: experiment with different values of pow, although
        pow=2 seems like a fine choice based on our limited
        experience."""
        return -((pl.absolute(pi - value) / sigma)**pow).sum()
    
    return vars()

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 build_moving5(days, avg):
    moving5 = array(zeros(len(days)-4), dtype = float)
    cday = 1
    moving5[0] = pylab.mean(avg[0:4])
    for a in avg[5:]:
	moving5[cday] = pylab.mean(avg[cday:cday+4])
	cday += 1
    return moving5

def perlin_covariance_corr(delta,N=1000000,bound=1):
    ts = bound*pl.rand(N)
    tds = ts+delta
    ps = [p(t) for t in ts]
    pds = [p(td) for td in tds]
    #cov = pl.mean([pp*pd for pp,pd in zip(ps,pds)])
    cov = pl.mean([(pp-pd)**2 for pp,pd in zip(ps,pds)])
    corr = pl.mean([pp*pd for pp,pd in zip(ps,pds)])
    return cov, corr

def int_f(a, fs=1.):
    """
    A fourier-based integrator.

    ===========
    Parameters:
    ===========
    a : *array* (1D)
        The array which should be integrated
    fs : *float*
        sampling time of the data

    ========
    Returns:
    ========
    y : *array* (1D)
        The integrated array

    """

    if False:
    # version with "mirrored" code
        xp = hstack([a, a[::-1]])
        int_fluc = int_f0(xp, float(fs))[:len(a)]
        baseline = mean(a) * arange(len(a)) / float(fs)
        return int_fluc + baseline - int_fluc[0]
    
    # old version
    baseline = mean(a) * arange(len(a)) / float(fs)
    int_fluc = int_f0(a, float(fs))
    return int_fluc + baseline - int_fluc[0]

    # old code - remove eventually (comment on 02/2014)
    # periodify
    if False:
        baseline = linspace(a[0], a[-1], len(a))
        a0 = a - baseline
        m = a0[-1] - a0[-2]
        b2 = linspace(0, -.5 * m, len(a))
        baseline -= b2
        a0 += b2
        a2 = hstack([a0, -1. * a0[1:][::-1]]) # "smooth" periodic signal  

        dbase = baseline[1] - baseline[0]
        t_vec = arange(len(a)) / float(fs)
        baseint = baseline[0] * t_vec + .5 * dbase * t_vec ** 2
        
        # define frequencies
        T = len(a2) / float(fs)
        freqs = 1. / T * arange(len(a2))
        freqs[len(freqs) // 2 + 1 :] -= float(fs)

        spec = fft.fft(a2)
        spec_i = zeros_like(spec, dtype=complex)
        spec_i[1:] = spec[1:] / (2j * pi* freqs[1:])
        res_int = fft.ifft(spec_i).real[:len(a0)] + baseint
        return res_int - res_int[0]

def xyamb(xytab,qu,xyout=''):

    mytb=taskinit.tbtool()

    if not isinstance(qu,tuple):
        raise Exception,'qu must be a tuple: (Q,U)'

    if xyout=='':
        xyout=xytab
    if xyout!=xytab:
        os.system('cp -r '+xytab+' '+xyout)

    QUexp=complex(qu[0],qu[1])
    print 'Expected QU = ',qu   # , '  (',pl.angle(QUexp)*180/pi,')'

    mytb.open(xyout,nomodify=False)

    QU=mytb.getkeyword('QU')['QU']
    P=pl.sqrt(QU[0,:]**2+QU[1,:]**2)

    nspw=P.shape[0]
    for ispw in range(nspw):
        st=mytb.query('SPECTRAL_WINDOW_ID=='+str(ispw))
        if (st.nrows()>0):
            q=QU[0,ispw]
            u=QU[1,ispw]
            qufound=complex(q,u)
            c=st.getcol('CPARAM')
            fl=st.getcol('FLAG')
            xyph0=pl.angle(pl.mean(c[0,:,:][pl.logical_not(fl[0,:,:])]),True)
            print 'Spw = '+str(ispw)+': Found QU = '+str(QU[:,ispw])  # +'   ('+str(pl.angle(qufound)*180/pi)+')'
            #if ( (abs(q)>0.0 and abs(qu[0])>0.0 and (q/qu[0])<0.0) or
            #     (abs(u)>0.0 and abs(qu[1])>0.0 and (u/qu[1])<0.0) ):
            if ( pl.absolute(pl.angle(qufound/QUexp)*180/pi)>90.0 ):
                c[0,:,:]*=-1.0
                xyph1=pl.angle(pl.mean(c[0,:,:][pl.logical_not(fl[0,:,:])]),True)
                st.putcol('CPARAM',c)
                QU[:,ispw]*=-1
                print '   ...CONVERTING X-Y phase from '+str(xyph0)+' to '+str(xyph1)+' deg'
            else:
                print '      ...KEEPING X-Y phase '+str(xyph0)+' deg'
            st.close()
    QUr={}
    QUr['QU']=QU
    mytb.putkeyword('QU',QUr)
    mytb.close()
    QUm=pl.mean(QU[:,P>0],1)
    QUe=pl.std(QU[:,P>0],1)
    Pm=pl.sqrt(QUm[0]**2+QUm[1]**2)
    Xm=0.5*atan2(QUm[1],QUm[0])*180/pi

    print 'Ambiguity resolved (spw mean): Q=',QUm[0],'U=',QUm[1],'(rms=',QUe[0],QUe[1],')','P=',Pm,'X=',Xm

    stokes=[1.0,QUm[0],QUm[1],0.0]
    print 'Returning the following Stokes vector: '+str(stokes)
    
    return stokes

def correctBias(AllData):
  # correct for difficulty and plot each subject %correct vs confidence
  corrmatrix, confmatrix = returnConfMatrix(AllData)
  Qs, subjects = py.shape(corrmatrix)
  copts = [1,2,3,4,5]
  datamat = np.array(py.zeros([len(copts), subjects]))
  print(datamat)
  fig = py.figure()
  ax15 = fig.add_subplot(111) 
  i = 0
 
  while i < subjects:
    c1, c2, c3, c4, c5 = [],[],[],[],[]
    # get confidences for each subject
    j = 0
    while j < Qs:
      # get confidences and correct for each question
      if confmatrix[j][i] == 1:
        c1.append(corrmatrix[j][i])
      elif confmatrix[j][i] == 2:
        c2.append(corrmatrix[j][i])
      elif confmatrix[j][i] == 3:
        c3.append(corrmatrix[j][i])
      elif confmatrix[j][i] == 4:
        c4.append(corrmatrix[j][i])
      elif confmatrix[j][i] == 5:
        c5.append(corrmatrix[j][i])
      else:
        print('bad num encountered')
        
      j += 1
    print('i is %d' %i)
    minconf = ([py.mean(c1), py.mean(c2), py.mean(c3), 
                   py.mean(c4), py.mean(c5)])
    pmin = 10
    for p in minconf:
      if p < pmin and p != 0 and math.isnan(p) is not True:
        pmin = p
    
    print(pmin)
    datamat[0][i] = py.mean(c1)/pmin
    datamat[1][i] = py.mean(c2)/pmin
    datamat[2][i] = py.mean(c3)/pmin
    datamat[3][i] = py.mean(c4)/pmin
    datamat[4][i] = py.mean(c5)/pmin
    # print(datamat)
    print( py.shape(datamat))
    print(len(datamat[:,i]))
    ax15.plot(range(1,6), datamat[:,i], alpha=0.4, linewidth=4)
    i += 1
  
  ax15.set_ylabel('Modified Correct')
  ax15.set_xlabel('Confidence')
  ax15.set_title('All responses')
  ax15.set_xticks(np.arange(1,6))
  ax15.set_xticklabels( [1, 2, 3, 4, 5] )
  ax15.set_xlim(0,6)

def nrms(data_fit, data_true):
    """
    Normalized root mean square error.
    """
    # root mean square error
    rms = pl.mean(pl.norm(data_fit - data_true, axis=0))

    # normalization factor is the max - min magnitude, or 2 times max dist from mean
    norm_factor = 2*pl.norm(data_true - pl.mean(data_true, axis=1), axis=0).max()
    return (norm_factor - rms)/norm_factor

 def ttest(X,Y):
     """
     Takes two lists of values, returns t value
     
     >>> ttest([2, 3, 7, 6, 10], [11,2,3,1,2])
     0.77459666924148329
     """
     if len(X) <= 1 or len(Y) <= 1: return 0.0
     return ((pylab.mean(X) - pylab.mean(Y))
             / stderr(X,Y))