Package gmisclib :: Module robust_multivariate
[frames] | no frames]

Source Code for Module gmisclib.robust_multivariate

  1   
  2  """Does a robust estimate of covariance. 
  3  It doesn't converge, unfortunately, for 
  4  vectors of dimension higher than 1. 
  5  """ 
  6   
  7  import Num 
  8  import math 
  9  import sys 
 10  import die 
 11   
12 -def cov_wt(veclist, weights):
13 nvec,vecdim = veclist.shape 14 o = Num.zeros((vecdim,vecdim), Num.Float) 15 for i in range(nvec): 16 Num.add(o, Num.outerproduct(veclist[i], veclist[i])*weights[i], o) 17 Num.divide(o, Num.sum(weights), o) 18 return o
19 20
21 -def integrate_guts(xlow, xhigh, fcn, c, n):
22 sum = 0.0 23 dx = (xhigh-xlow)/float(n) 24 for i in range(n): 25 x = xlow + (i+0.5)*dx 26 sum += fcn(x, c) 27 return sum*dx
28 29
30 -def integrate(xlow, xhigh, fcn, c=None, nstart=5, tol=1e-4):
31 """A generally useful function, does the integral of 32 fcn(x, c) from xlow to xhigh. 33 """ 34 passes = 0 35 while passes < 10: 36 i = integrate_guts(xlow, xhigh, fcn, c, nstart) 37 tmp = integrate_guts(xlow, xhigh, fcn, c, nstart*3) 38 if abs(tmp-i) < tol: 39 return tmp + (tmp-i)/9.0 40 nstart *= 3 41 i = tmp 42 passes += 1 43 raise RuntimeError, "Integral doesn't converge"
44 45
46 -def wtfunc(q, E):
47 Qcut = 2.0 48 return (Qcut/(Qcut + q))**E
49 # return Num.minimum(1.0, q**(-E)) 50 51
52 -def wtdim(q, Edim):
53 E,vecdim,qxtra = Edim 54 q2 = q**2 55 return math.exp(-q2/2.0)*wtfunc(q2,E) * q**qxtra
56 57
58 -def covariance(veclist, emax=1.0):
59 """Veclist[i] is a vector. 60 """ 61 echg = 0.25 62 E = echg * emax 63 nvec, vecdim = veclist.shape 64 e0 = integrate(0, 3.0+0.5*vecdim, wtdim, c=(0.0,vecdim, 0), 65 nstart = vecdim*2) 66 e0qq = integrate(0, 3.0+0.5*vecdim, wtdim, c=(0.0,vecdim, 2), 67 nstart = vecdim*2) 68 # Qtol = SMALL * NDIM * WIDTH_OF_CHISQ(ndim)_DISTRIBUTION 69 qtol = 0.001*nvec*math.sqrt(vecdim) 70 weights = Num.ones((len(veclist),), Num.Float) 71 efac = 1.0 72 qlast = None 73 ipass = 0 74 while ipass < 1000: 75 c = cov_wt(veclist, weights) / efac 76 assert c.shape == (vecdim, vecdim) 77 print 'c.shape=', c.shape 78 print 'c=', c 79 eval, evec = Num.LA.eigh(c) 80 print 'evals=', eval 81 if not Num.alltrue(eval >= 0.0): 82 weights = Num.sqrt(weights) 83 die.warn('Negative eigenvalue: slowing down.') 84 continue 85 ichalf = evec/Num.sqrt(eval) 86 # print 'ERR=', Num.dot(Num.transpose(ichalf), Num.dot(c, ichalf)) 87 # print 'ichalf=', ichalf 88 vecxf2 = Num.dot(veclist, ichalf)**2 89 print 'vecxf2=', vecxf2 90 assert vecxf2.shape == (nvec, vecdim) 91 weighttmp = wtfunc(vecxf2, E) 92 print 'weighttmp=', weighttmp 93 weights = Num.prod(weighttmp, axis=1) 94 assert weights.shape == (nvec,) 95 q = Num.sum(vecxf2, axis=1) 96 assert q.shape == (nvec,) 97 print 'q=', q 98 assert Num.alltrue(q >= 0.0) 99 print 'weight=', weights 100 if qlast is not None: 101 qchg = Num.sum(Num.absolute(q-qlast)) 102 echg = 1.0/math.sqrt(qchg/qtol + 1.0) 103 print 'qchg=', qchg, 'echg=', echg, 'E=', E 104 if qchg < qtol and abs(E-emax)<0.01: 105 return (c, eval, evec, q) 106 qlast = q 107 e = integrate(0, 3.0+0.5*vecdim, wtdim, c=(E,vecdim,0), 108 nstart = vecdim*2) 109 eqq = integrate(0, 3.0+0.5*vecdim, wtdim, c=(E,vecdim,2), 110 nstart = vecdim*2) 111 efac = (eqq/e0qq)/(e/e0) 112 print 'efac=', efac 113 E = echg*emax + (1.0-echg)*E 114 ipass += 1 115 sys.stdout.flush() 116 raise RuntimeError, "Covariance does not converge."
117 118 119 if __name__ == '__main__': 120 Nvec = 2000 121 ND = 10 122 veclist = Num.RA.normal(0.0, 1.0, (Nvec,ND)) 123 covariance(veclist) 124