Source Code for Module classifiers.multivariance

  1  from gmisclib import Num 
  2  from gmisclib import mcmc 
  3  from gmisclib import die 
  4  from gmisclib import mcmc_helper 
  5  import gmisclib.multivariance_classes as MC 
  6  import gmisclib.multivariance_q as MQ 
  7  import gmisclib.multivariance_mm as MM 
  8   
  9   
 10   
 11  ERGCOVER = 3.0 
 12  HUGE = 1e30 
 13                                  # -1 for a flat prior in inverse sigma, 
 14                                  # 0 for somewhere in between 
 15   
 16   
 17   
 18   
 19   
 20 -def _logp(prms, (x, modelchoice)): 
 21          lp = 0.0 
 22          p = modelchoice.unpack(prms) 
 23   
 24          try: 
 25                  for x1 in x: 
 26                          lp += p.logp(x1) 
 27          except MC.QuadraticNotNormalizable: 
 28                  return -HUGE 
 29   
 30          return lp 
 31   
 32   
 33   
 34   
 35   
 36 -def meanvar(dataset, N, modelchoice=MQ.quadratic): 
 37          """Given a dataset, this produces a bunch of MCMC estimates 
 38          of plausible means and inverse covariance matrices for the dataset. 
 39          Modelchoice is a particular size model, e.g. an instance of quadratic. 
 40          The output is in the form [ (mean, inv_covar), ... ].""" 
 41          assert N >= 0 
 42          assert len(dataset)>0 
 43          # print 'N=', N 
 44          particular_model = modelchoice(dataset=dataset) 
 45          startpack, varpack = particular_model.start(dataset) 
 46          # print 'SPACK:', startpack 
 47          # print 'VPACK:', varpack 
 48          x = mcmc.bootstepper(_logp, startpack, varpack, 
 49                                  c=(dataset,particular_model)) 
 50          # print 'X!' 
 51   
 52          n = particular_model.ndim() 
 53          # print 'ndim=', n 
 54          mcmch = mcmc_helper.stepper(x) 
 55          mcmch.run_to_bottom() 
 56          # print 'bottom', N 
 57          o = [] 
 58          for i in range(N): 
 59                  mcmch.run_to_ergodic(ERGCOVER/float(N)) 
 60                  o.append( particular_model.unpack(x.prms()) ) 
 61          return o 
 62   
 63   
 64 -def test_quadratic(): 
 65          d = [] 
 66          for l in sys.stdin: 
 67                  if l.startswith('#'): 
 68                          continue 
 69                  if l.strip() == '': 
 70                          continue 
 71                  a = Num.array( [float(s) for s in l.split()], Num.Float) 
 72                  d.append(a) 
 73          for x in meanvar(d, 200, modelchoice=MQ.quadratic): 
 74                  print 'Mu:', x.mu 
 75                  print 'Inverse:', x.invsigma 
 76                  print 
 77   
 78   
 79 -def test_diag_quadratic(): 
 80          d = [] 
 81          for l in sys.stdin: 
 82                  if l.startswith('#'): 
 83                          continue 
 84                  if l.strip() == '': 
 85                          continue 
 86                  a = Num.array( [float(s) for s in l.split()], Num.Float) 
 87                  d.append(a) 
 88          for x in meanvar(d, 200, modelchoice=MQ.diag_quadratic): 
 89                  print 'Mu:', x.mu 
 90                  print 'Inverse:', x.invsigma 
 91                  print 
 92   
 93 -def test_multimu(): 
 94          d = [] 
 95          for l in sys.stdin: 
 96                  if l.startswith('#'): 
 97                          continue 
 98                  if l.strip() == '': 
 99                          continue 
100                  a = l.split() 
101                  af = Num.array([float(s) for s in a[1:]], Num.Float) 
102                  d.append(MM.datum_c(vector=af, classid=a[0])) 
103          for x in meanvar(d, 200, modelchoice=MM.multi_mu): 
104                  print 'Inverse:', x.invsigma 
105                  print 'Mu:', x.mu 
106                  print 
107   
108   
109 -def test_multimu_diag(): 
110          d = [] 
111          for l in sys.stdin: 
112                  if l.startswith('#'): 
113                          continue 
114                  if l.strip() == '': 
115                          continue 
116                  a = l.split() 
117                  af = Num.array([float(s) for s in a[1:]], Num.Float) 
118                  d.append(MM.datum_c(vector=af, classid=a[0])) 
119          for x in meanvar(d, 200, modelchoice=MM.multi_mu_diag): 
120                  print 'Inverse:', x.invsigma 
121                  print 'Mu:', x.mu 
122                  print 
123   
124   
125  FromName = {'quadratic' : MQ.quadratic, 'diag_quadratic': MQ.diag_quadratic, 
126                  'multi_mu': MM.multi_mu, 'multi_mu_diag': MM.multi_mu_diag 
127          } 
128   
129  if __name__ == '__main__': 
130          import sys 
131   
132          try: 
133                  import psyco 
134                  psyco.full(memory=10000) 
135          except ImportError: 
136                  pass 
137   
138          if len(sys.argv) <= 1: 
139                  print 'Usage: multivariance model_name ...' 
140                  print __doc__ 
141                  die.exit(0) 
142          if sys.argv[1] == 'quadratic': 
143                  sys.argv = [sys.argv[0]] + sys.argv[2:] 
144                  test_quadratic() 
145          elif sys.argv[1] == 'diag_quadratic': 
146                  sys.argv = [sys.argv[0]] + sys.argv[2:] 
147                  test_diag_quadratic() 
148          elif sys.argv[1] == 'multi_mu': 
149                  sys.argv = [sys.argv[0]] + sys.argv[2:] 
150                  test_multimu() 
151          elif sys.argv[1] == 'multi_mu_diag': 
152                  sys.argv = [sys.argv[0]] + sys.argv[2:] 
153                  test_multimu_diag() 
154