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
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
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
47 Qcut = 2.0
48 return (Qcut/(Qcut + q))**E
49
50
51
53 E,vecdim,qxtra = Edim
54 q2 = q**2
55 return math.exp(-q2/2.0)*wtfunc(q2,E) * q**qxtra
56
57
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
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
87
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