/home/mark/model/software/ScrumPy/ScrumPy/Data/Stats.py

import math,random

#import sys
#sys.path.extend(["../ThirdParty", "../Util"])    ####### testing only
#import Special,  Sci

from ScrumPy.Util import  Sci
from ScrumPy.ThirdParty import Special

Pi2 = 2*math.pi
RootPi2 = math.sqrt(Pi2)
Root2 = math.sqrt(2.0)
NaN=float("NaN")

def NormPdf(x,mu=0,sigma=1):
    frac = (x-mu)/sigma
    Exp = math.exp(-0.5*frac*frac)
    recip = 1/(sigma*RootPi2)
    return recip * Exp


def NormCdf(x,mu=0,sigma=1):

    erfarg = (x-mu)/(sigma*Root2)
    return (1+Special.Erf(erfarg))/2



def Binv(p,q):
    """inverse beta function"""
    # from, and tested against, gnuplot stat.inc
    lgamma = Special.Lgamma
    return math.exp(lgamma(p+q)-lgamma(p)-lgamma(q))



def t_pdf(x,nu):
    """ student's t pdf for nu d.f. """
    # from, and tested against, gnuplot stat.inc
    #t(x,nu)=nu<0||!isint(nu)?1/0:Binv(0.5*nu,0.5)/sqrt(nu)*(1+real(x*x)/nu)**(-0.5*(nu+1.0))

    return Binv(0.5*nu,0.5)/math.sqrt(nu)*(1+(x*x)/nu)**(-0.5*(nu+1.0))



def t_cdf(x,nu):
    """ student's t cdf for nu d.f. """
    #  from, and tested against, gnuplot stat.inc
    # ct(x,nu)=nu<0||!isint(nu)?1/0:0.5+0.5*sgn(x)*(1-ibeta(0.5*nu,0.5,nu/(nu+x*x)))
    ibeta = Special.betai

    if x<0:
        sgn5 = -0.5
    else:
        sgn5 = 0.5

    return 0.5+sgn5*(1-ibeta(0.5*nu,0.5,nu/(nu+x*x)))


def inv_t_cdf(c,nu, tol=1e-6, lim=30):

    # tol and lim good for at least .5025 <= c <=1-1e-7, n >= 1
    # (opteron h/w)

    def itcdf(lo,hi,n=0):

        mid = (lo+hi)/2
        cur = t_cdf(mid,nu)
        err = 1 - lo/mid

        if err <= tol:
            return mid

        if n>lim:
            print "!! itcdf: only achieved tol=", err, "requested", tol, "!!"
            return mid

        if cur > c:
            return itcdf(lo,mid,n+1)

        return itcdf(mid,hi,n+1)

    inihi = 5.
    while t_cdf(inihi,nu)<c:
        inihi *=2


    return itcdf(0.0,inihi)


def Conflims(conf, s, df):
    # s is sample sd

    # see Brown & Rothery (1ed) 5.4.3
    c = 1 - (1-conf)/2
    tval =  inv_t_cdf(c,df-1) # c/w B&R table S2
    return tval *s/math.sqrt(df)










def NormList(mu,sigma,size):

    return map(lambda x: random.normalvariate(mu,sigma), [0]*size)


def SumProds(x,y):
    """ pre: len(x) == len(y), items are numeric
      post: sum(x[i] * y[i]) : i = 0..len(x)-1 """
    return sum(map(lambda a,b: a*b, x,y))

def SqSum(x):
    """ pre: sum(x)
       post: sum(x)**2 """

    s = sum(x)
    return s*s

def SumSq(x):
    return sum(map(lambda a:a*a,x))

def SumSqDiff(vals, const):

        return(SumSq(map(lambda x: x-const, vals)))


def Adj_p(plist):
    rv = plist[:]
    rv.sort()
    mul = len(rv)
    for i in range(len(plist)):
        rv[i] *= mul
        mul -= 1
    return rv

def sortidx(x):
    """
    faster equivalent of stats.lshellsort(inlist)
    """
    n = len(x)
    xs,x = x[:],x[:]
    xs.sort()
    rv = [0] *n
    for i in range(n):
            rv[i] = x.index(xs[i])
            x[rv[i]] = None
    return xs,rv


def Rank(inlist):
    """
Ranks the data in inlist, dealing with ties appropritely.  Assumes
a 1D inlist.  Adapted from Gary Perlman's |Stat ranksort.

Usage:   lrankdata(inlist)
Returns: a list of length equal to inlist, containing rank scores
"""
    n = len(inlist)
    svec, ivec = sortidx(inlist)
    sumranks = 0
    dupcount = 0
    newlist = [0]*n
    for i in range(n):
        sumranks +=  i
        dupcount +=  1
        if i==n-1 or svec[i] <> svec[i+1]:
            averank = sumranks / float(dupcount )+ 1
            for j in range(i-dupcount+1,i+1):
                newlist[ivec[j]] = averank
            sumranks = 0
            dupcount = 0
    return newlist






def Mean(a):
    return sum(a)/len(a)

def Var(a,abar=None):
    if abar ==None:
        abar = Mean(a)

    return sum( map(lambda x:(x-abar)**2,a))/(len(a)-1)


def StdDev(a, abar = None):

    if abar == None:
        abar = Mean(a)

    return math.sqrt(Var(a,abar))


def StdErr(a):

        return StdDev(a)/math.sqrt(len(a))


def Skewness(a,mu=None,sigma=None): # NRiC 13.1.5

    lena = len(a)
    if mu==None:
        mu = Mean(a)
    if sigma==None:
        sigma = StdDev(a,mu)

    return sum(  map(lambda x: ((x - mu)/sigma)**3,a))/lena


def Kurtosis(a):# NRiC 13.1.6
    lena = len(a)
    mean = Mean(a)
    sigma = StdDev(a)

    return (sum(  map(lambda x: ((x - mean)/sigma)**4,a))/lena)-3



def Lower_n_ile(vals, n,UseCopy=True):
    """ pre: len(self[Name]) > 0, 0 < n <= 1
       post: self.Lower_n_ile(self, Name, n) == value of self[Name] which is the limit of the
             lower nth portion of self[Name] """

    if UseCopy:
        vals = vals[:]
    vals.sort()
    index = int(len(vals)*n)
    return vals[index]


def Upper_n_ile(vals, n,UseCopy=True):
    """ pre: len(self[Name]) > 0, 0 < n <= 1
       post: self.Lower_n_ile(self, Name, n) == value of self[Name] which is the limit of the
             upper nth portion of self[Name] """

    return Lower_n_ile(vals, 1.0-n,UseCopy)

def Median(vals, UseCopy=True):
    """ pre: len(self[Name]) > 0
        post: self.Median(Name) == Median value of self[Name] """

    return Lower_n_ile(vals, 0.5,UseCopy)



def FTest(vals1, vals2): # NRiC 13.4
   """ F and p(F) for unequal variance
   """

   var1 = Var(vals1)
   var2 = Var(vals2)
   df1 = len(vals1)-1
   df2 = len(vals2)-1

   f = var1/var2
   if f<1:
           f  = 1/f
           df1, df2 = df2, df1

   p = 2* Special.betai(0.5*df2, 0.5*df1, df2/(f*df1+df2))

   if p>1:
           p = 2-p

   return f, p


def TTest(vals1, vals2):
        """ Student's (two tailed) t and p(t) fo different means
        """

        # NRiC 13.4.3/4
        # assume unequal var. Equivalent if vars are equal anyway

        v1bar = Mean(vals1)
        v2bar = Mean(vals2)
        v1var = Var(vals1)
        v2var = Var(vals2)
        n1 = len(vals1)
        n2 = len(vals2)

        t = (v1bar-v2bar)/(
             math.sqrt(
                v1var/n1 + v2var/n2
                )
            )

        df = (v1var/n1 + v2var/n2)**2/(
                ((v1var/n1)**2/(n1-1)) +
                ((v2var/n2)**2/(n2-1))
            )


        return t, Prob_of_t(t, df)


def Prob_of_t(t, df):

        df = float(df)

        return Special.betai(0.5*df, 0.5, df/(df+t*t))





def CoVar(a,b,abar=None, bbar=None):

    if abar == None: abar = Mean(a)
    if bbar == None: bbar = Mean(b)

    return   sum(map(lambda x,y:(x-abar)*(y-bbar),a,b))/(len(a)-1)


def Pearson_r(x,y):
    """
    pre: len(x) == len(y),  items are numeric type
    post: Pearsons r
    """

    xmean = Mean(x)
    ymean = Mean(y)
    xsubmean = map(lambda a: a-xmean,x)
    ysubmean = map(lambda a: a-ymean, y)
    sq = lambda n: n*n

    r_num = SumProds(xsubmean, ysubmean)
    r_den = math.sqrt(sum(map(sq, xsubmean))) * math.sqrt(sum(map(sq, ysubmean)))

    if r_den == 0:
        return NaN

    return r_num / r_den



def Spearman_r(x, y, rankxy=None,):

    """
Calculates a Spearman rank-order correlation coefficient.  Taken
from Heiman's Basic Statistics for the Behav. Sci (1st), p.192.

Usage:   lspearmanr(x,y)      where x and y are equal-length lists
Returns: Spearman's r, two-tailed p-value
"""

    n = len(x)
    if rankxy == None:
            rankx = Rank(x)
            ranky = Rank(y)
    else:
        rankx, ranky = rankxy

    diffs = map(lambda a,b: a-b, rankx,ranky)
    dsq = SumSq(diffs)
    return 1 - 6*dsq / float(n*(n**2-1))

def Prob_of_r(r, n):
    df = n-2
    t = r * math.sqrt(df / ( (r+1.0)*(1.0-r)))
    return Special.betai(0.5*df, 0.5,  df/(df+t*t) )





def Prob_of_chisq(chisq,n):
    """ probability of getting <chisq for correct model with n deg freedom.
        NRiC 6.2.15 """

    return Special.Gamma_p(float(n),chisq)




"""
Calculates a one-way chi square for list of observed frequencies and returns
the result.  If no expected frequencies are given, the total N is assumed to
be equally distributed across all groups.

Usage:   lchisquare(Samp, f_exp=None)   f_obs = list of observed cell freq.
Returns: chisquare-statistic, associated p-value
"""

def ChiSquare(Samp, Ref, NormRef=True):

    if NormRef:
        coeff = float(sum(Ref))/sum(Samp)
        Ref = map(lambda x:x/coeff,Ref)

    return sum(map(lambda s,r:((s-r)**2)/r,Samp,Ref))


def LinFit(Xseq, Yseq, Resids=False):
    """
    pre: Xseq, Yseq equal length sequence of numbers

    Calculate best m,c for y = mx+c

    post:
       Resids=False : return (m,c)
               else : return (m,c, [resids])
    """

    m,c = Sci.polyfit(Xseq,Yseq)
    if Resids:
        best = map(lambda x:m*x+c, Xseq)
        resids = map(lambda x,y:x-y, best, Yseq)
        rv = m,c,resids
    else:
        rv = m,c

    return rv



##def ClustTest(List, CritFun, binsize, nulchecks):
##
##    hits = map(CritFun, List)
##    print hits.count(True),
##    hhist = BoolList2His(hits,binsize)
##    Chi_p = stats.chisquare(*Lis(hhist))[1]
##    nulps = []
##    for n in range(nulchecks):
##        random.shuffle(hits)
##        nhist = BooList2His(hits,binsize)
##        nulps.append(stats.chisquare(HisToLis(nhist)))
##    return Chi_p,nulps



def ColMeansMtx(mtx):

    return map(lambda c: Mean(mtx.GetCol(c)),mtx.cnames)


def CoVarMtx(mtx):

    lenr = len(mtx.rnames)
    coeff = math.sqrt(lenr-1.)
    mtx = mtx.Copy()
    negmeans = map(lambda x:-x,ColMeansMtx(mtx))

    for r in range(lenr):
        mtx.AddRow(r, negmeans)
        mtx.DivRow(r,k=coeff)

    return mtx.Copy(tx=1).Mul(mtx)


def CorMtx(mtx):
    raise  "use pearson instead"

def Rank_mtx(mtx):
    rv = mtx.Copy()
    for c in mtx.cnames:
        rv.SetCol(c, Rank(mtx.GetCol(c)))
    return rv

def Spearman_rmtx(mtx):
    rmtx = Rank_mtx(mtx) # make a cache of the ranks, save recomputing
    rmtx.Transpose()

    def spear(m,a,b):
        return Spearman_r(m[a],m[b],rankxy=(rmtx[a],rmtx[b]))

    mtx.Transpose()
    rv = mtx.RowDiffMtx2(spear)
    for n in range(len(rv)):
        rv[n,n] = 1.0
    mtx.Transpose()
    return rv

def Pearson_rmtx(mtx):
    mtx.Transpose()
    rv = mtx.RowDiffMtx(Pearson_r)
    mtx.Transpose()
    for n in range(len(rv)):
        rv[n,n] = 1.0
    return rv

def Prob_rmtx(rmtx,n):
    """
    pre: rmtx is  a matrix of r values, n>2 the number of samples used to determine it
    post: matrix of corresponding p values
    """
    rv = rmtx.Copy()
    n = len(rmtx)
    for r1 in range(n):
        for r2 in range(r1+1,n):
            rv[r1,r2] = rv[r2,r1] = Prob_of_r(rmtx[r1,r2],n)

    for r in range(n):
        rv[r,r] = 0.0

    return rv