Source Code for Module rdkit.ML.Data.Stats

  1  ## Automatically adapted for numpy.oldnumeric Jun 27, 2008 by -c 
  2   
  3  # $Id: Stats.py 997 2009-02-25 06:12:43Z glandrum $ 
  4  # 
  5  #  Copyright (C) 2001-2008  greg Landrum and Rational Discovery LLC 
  6  #  All Rights Reserved 
  7  # 
  8  """ various statistical operations on data 
  9   
 10  """ 
 11  import numpy 
 12   
 13 -def StandardizeMatrix(mat): 
 14    """ 
 15   
 16    This is the standard *subtract off the average and divide by the deviation* 
 17    standardization function. 
 18   
 19     **Arguments** 
 20   
 21       - mat: a numpy array 
 22   
 23     **Notes** 
 24   
 25       - in addition to being returned, _mat_ is modified in place, so **beware** 
 26   
 27    """ 
 28    nObjs = len(mat) 
 29    avgs = sum(mat,0)/float(nObjs) 
 30    mat -= avgs 
 31    devs =sqrt(sum(mat*mat,0)/(float(nObjs-1))) 
 32    try: 
 33      newMat = mat/devs 
 34    except OverflowError: 
 35      newMat = numpy.zeros(mat.shape,'d') 
 36      for i in range(mat.shape[1]): 
 37        if devs[i] != 0.0: 
 38          newMat[:,i] = mat[:,i]/devs[i] 
 39    return newMat 
 40   
 41 -def FormCovarianceMatrix(mat): 
 42    """ form and return the covariance matrix 
 43   
 44    """ 
 45    nPts = mat.shape[0] 
 46    sumVect = sum(mat) 
 47    sumVect /= float(nPts) 
 48    for row in mat: 
 49      row -= sumVect 
 50    return numpy.dot(numpy.transpose(mat),mat)/(nPts-1) 
 51   
 52 -def FormCorrelationMatrix(mat): 
 53    """ form and return the covariance matrix 
 54   
 55    """ 
 56    nVars = len(mat[0]) 
 57    N = len(mat) 
 58     
 59    res = numpy.zeros((nVars,nVars),'d') 
 60    for i in range(nVars): 
 61      x = mat[:,i] 
 62      sumX = sum(x) 
 63      sumX2 = sum(x*x) 
 64      for j in range(i,nVars): 
 65        y = mat[:,j] 
 66        sumY = sum(y) 
 67        sumY2 = sum(y*y) 
 68        numerator = N*sum(x*y) - sumX*sumY 
 69        denom = numpy.sqrt((N*sumX2-sumX**2)*(N*sumY2-sumY**2)) 
 70        if denom != 0.0: 
 71          res[i,j] = numerator/denom 
 72          res[j,i] = numerator/denom 
 73        else: 
 74          res[i,j] = 0 
 75          res[j,i] = 0 
 76    return res 
 77 -def PrincipalComponents(mat,reverseOrder=1): 
 78    """ do a principal components analysis 
 79   
 80    """ 
 81    import numpy.oldnumeric.linear_algebra as LinearAlgebra 
 82    covMat = FormCorrelationMatrix(mat) 
 83   
 84    eigenVals,eigenVects = LinearAlgebra.eigenvectors(covMat) 
 85    try: 
 86      eigenVals = eigenVals.real 
 87    except: 
 88      pass 
 89    try: 
 90      eigenVects = eigenVects.real 
 91    except: 
 92      pass 
 93   
 94    # and now sort: 
 95    ptOrder = numpy.argsort(eigenVals).tolist() 
 96    if reverseOrder: 
 97      ptOrder.reverse() 
 98    eigenVals = numpy.array([eigenVals[x] for x in ptOrder]) 
 99    eigenVects = numpy.array([eigenVects[x] for x in ptOrder]) 
100    return eigenVals,eigenVects 
101   
102 -def TransformPoints(tFormMat,pts): 
103    """ transforms a set of points using tFormMat 
104   
105      **Arguments** 
106   
107        - tFormMat: a numpy array 
108   
109        - pts: a list of numpy arrays (or a 2D array) 
110   
111      **Returns** 
112   
113        a list of numpy arrays  
114   
115    """ 
116    pts = numpy.array(pts) 
117    nPts = len(pts) 
118    avgP = sum(pts)/nPts 
119    pts = pts - avgP 
120    res = [None]*nPts 
121    for i in range(nPts): 
122      res[i] = numpy.dot(tFormMat,pts[i]) 
123   
124    return res   
125   
126 -def MeanAndDev(vect,sampleSD=1): 
127    """ returns the mean and standard deviation of a vector """ 
128    vect = numpy.array(vect,'d') 
129    n = vect.shape[0] 
130    if n <= 0: 
131      return 0.,0. 
132    mean = sum(vect)/n 
133    v = vect-mean 
134    if n > 1: 
135      if sampleSD: 
136        dev = numpy.sqrt(sum(v*v)/(n-1)) 
137      else: 
138        dev = numpy.sqrt(sum(v*v)/(n)) 
139   
140    else: 
141      dev = 0 
142    return mean,dev 
143   
144 -def R2(orig,residSum): 
145    """ returns the R2 value for a set of predictions """ 
146   
147   
148    # FIX: this just is not right  
149    # 
150    #  A correct formulation of this (from Excel) for 2 variables is: 
151    #   r2 = [n*(Sxy) - (Sx)(Sy)]^2 / ([n*(Sx2) - (Sx)^2]*[n*(Sy2) - (Sy)^2]) 
152    # 
153    # 
154   
155    vect = numpy.array(orig) 
156    n = vect.shape[0] 
157    if n <= 0: 
158      return 0.,0. 
159    oMean = sum(vect)/n 
160    v = vect-oMean 
161    oVar = sum(v*v) 
162    return 1. - residSum/oVar 
163   
164   
165  # One Tail   0.10    0.05    0.025   0.01    0.005   0.001  0.0005 
166  tConfs = {80:1,90:2,95:3,98:4,99:5,99.8:6,99.9:7} 
167  tTable=[ 
168  [  1,    3.078,   6.314,   12.71,   31.82,   63.66,  318.30,     637], 
169  [  2,    1.886,   2.920,   4.303,   6.965,   9.925,  22.330,    31.6], 
170  [  3,    1.638,   2.353,   3.182,   4.541,   5.841,  10.210,   12.92], 
171  [  4,    1.533,   2.132,   2.776,   3.747,   4.604,   7.173,   8.610], 
172  [  5,    1.476,   2.015,   2.571,   3.365,   4.032,   5.893,   6.869], 
173  [  6,    1.440,   1.943,   2.447,   3.143,   3.707,   5.208,   5.959], 
174  [  7,    1.415,   1.895,   2.365,   2.998,   3.499,   4.785,   5.408], 
175  [  8,    1.397,   1.860,   2.306,   2.896,   3.355,   4.501,   5.041], 
176  [  9,    1.383,   1.833,   2.262,   2.821,   3.250,   4.297,   4.781], 
177  [ 10,    1.372,   1.812,   2.228,   2.764,   3.169,   4.144,   4.587], 
178  [ 11,    1.363,   1.796,   2.201,   2.718,   3.106,   4.025,   4.437], 
179  [ 12,    1.356,   1.782,   2.179,   2.681,   3.055,   3.930,   4.318], 
180  [ 13,    1.350,   1.771,   2.160,   2.650,   3.012,   3.852,   4.221], 
181  [ 14,    1.345,   1.761,   2.145,   2.624,   2.977,   3.787,   4.140], 
182  [ 15,    1.341,   1.753,   2.131,   2.602,   2.947,   3.733,   4.073], 
183  [ 16,    1.337,   1.746,   2.120,   2.583,   2.921,   3.686,   4.015], 
184  [ 17,    1.333,   1.740,   2.110,   2.567,   2.898,   3.646,   3.965], 
185  [ 18,    1.330,   1.734,   2.101,   2.552,   2.878,   3.610,   3.922], 
186  [ 19,    1.328,   1.729,   2.093,   2.539,   2.861,   3.579,   3.883], 
187  [ 20,    1.325,   1.725,   2.086,   2.528,   2.845,   3.552,   3.850], 
188  [ 21,    1.323,   1.721,   2.080,   2.518,   2.831,   3.527,   3.819], 
189  [ 22,    1.321,   1.717,   2.074,   2.508,   2.819,   3.505,   3.792], 
190  [ 23,    1.319,   1.714,   2.069,   2.500,   2.807,   3.485,   3.768], 
191  [ 24,    1.318,   1.711,   2.064,   2.492,   2.797,   3.467,   3.745], 
192  [ 25,    1.316,   1.708,   2.060,   2.485,   2.787,   3.450,   3.725], 
193  [ 26,    1.315,   1.706,   2.056,   2.479,   2.779,   3.435,   3.707], 
194  [ 27,    1.314,   1.703,   2.052,   2.473,   2.771,   3.421,   3.690], 
195  [ 28,    1.313,   1.701,   2.048,   2.467,   2.763,   3.408,   3.674], 
196  [ 29,    1.311,   1.699,   2.045,   2.462,   2.756,   3.396,   3.659], 
197  [ 30,    1.310,   1.697,   2.042,   2.457,   2.750,   3.385,   3.646], 
198  [ 32,    1.309,   1.694,   2.037,   2.449,   2.738,   3.365,   3.622], 
199  [ 34,    1.307,   1.691,   2.032,   2.441,   2.728,   3.348,   3.601], 
200  [ 36,    1.306,   1.688,   2.028,   2.434,   2.719,   3.333,   3.582], 
201  [ 38,    1.304,   1.686,   2.024,   2.429,   2.712,   3.319,   3.566], 
202  [ 40,    1.303,   1.684,   2.021,   2.423,   2.704,   3.307,   3.551], 
203  [ 42,    1.302,   1.682,   2.018,   2.418,   2.698,   3.296,   3.538], 
204  [ 44,    1.301,   1.680,   2.015,   2.414,   2.692,   3.286,   3.526], 
205  [ 46,    1.300,   1.679,   2.013,   2.410,   2.687,   3.277,   3.515], 
206  [ 48,    1.299,   1.677,   2.011,   2.407,   2.682,   3.269,   3.505], 
207  [ 50,    1.299,   1.676,   2.009,   2.403,   2.678,   3.261,   3.496], 
208  [ 55,    1.297,   1.673,   2.004,   2.396,   2.668,   3.245,   3.476], 
209  [ 60,    1.296,   1.671,   2.000,   2.390,   2.660,   3.232,   3.460], 
210  [ 65,    1.295,   1.669,   1.997,   2.385,   2.654,   3.220,   3.447], 
211  [ 70,    1.294,   1.667,   1.994,   2.381,   2.648,   3.211,   3.435], 
212  [ 80,    1.292,   1.664,   1.990,   2.374,   2.639,   3.195,   3.416], 
213  [100,    1.290,   1.660,   1.984,   2.364,   2.626,   3.174,   3.390], 
214  [150,    1.287,   1.655,   1.976,   2.351,   2.609,   3.145,   3.357], 
215  [200,    1.286,   1.653,   1.972,   2.345,   2.601,   3.131,   3.340] 
216  ] 
217 -def GetConfidenceInterval(sd,n,level=95): 
218    col = tConfs[level] 
219    dofs = n-1 
220    sem = sd/sqrt(n) 
221    idx = 0 
222    while idx<len(tTable) and tTable[idx][0]<dofs: 
223      idx+=1 
224    if idx<len(tTable): 
225      t = tTable[idx][col] 
226    else: 
227      t = tTable[-1][col] 
228    return t*sem 
229