[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):
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
<!--
expandto(location.href);
// -->

```

 Generated by Epydoc 3.0.1 on Thu Sep 22 04:25:15 2011 http://epydoc.sourceforge.net