Source Code for Module mvpa.clfs.plr

  1  # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- 
  2  # vi: set ft=python sts=4 ts=4 sw=4 et: 
  3  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 
  4  # 
  5  #   See COPYING file distributed along with the PyMVPA package for the 
  6  #   copyright and license terms. 
  7  # 
  8  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 
  9  """Penalized logistic regression classifier.""" 
 10   
 11  __docformat__ = 'restructuredtext' 
 12   
 13   
 14  import numpy as N 
 15   
 16  from mvpa.misc.exceptions import ConvergenceError 
 17  from mvpa.clfs.base import Classifier, FailedToTrainError 
 18   
 19  if __debug__: 
 20      from mvpa.base import debug 
 21   
 22   
 23 -class PLR(Classifier): 
 24      """Penalized logistic regression `Classifier`. 
 25      """ 
 26   
 27      _clf_internals = [ 'plr', 'binary', 'linear' ] 
 28   
 29 -    def __init__(self, lm=1, criterion=1, reduced=0.0, maxiter=20, **kwargs): 
 30          """ 
 31          Initialize a penalized logistic regression analysis 
 32   
 33          :Parameters: 
 34            lm : int 
 35              the penalty term lambda. 
 36            criterion : int 
 37              the criterion applied to judge convergence. 
 38            reduced : float 
 39              if not 0, the rank of the data is reduced before 
 40              performing the calculations. In that case, reduce is taken 
 41              as the fraction of the first singular value, at which a 
 42              dimension is not considered significant anymore. A 
 43              reasonable criterion is reduced=0.01 
 44            maxiter : int 
 45              maximum number of iterations. If no convergence occurs 
 46              after this number of iterations, an exception is raised. 
 47   
 48          """ 
 49          # init base class first 
 50          Classifier.__init__(self, **kwargs) 
 51   
 52          self.__lm   = lm 
 53          self.__criterion = criterion 
 54          self.__reduced = reduced 
 55          self.__maxiter = maxiter 
 56   
 57   
 58 -    def __repr__(self): 
 59          """String summary over the object 
 60          """ 
 61          return """PLR(lm=%f, criterion=%d, reduced=%s, maxiter=%d, enable_states=%s)""" % \ 
 62                 (self.__lm, self.__criterion, self.__reduced, self.__maxiter, 
 63                  str(self.states.enabled)) 
 64   
 65   
 66 -    def _train(self, data): 
 67          """Train the classifier using `data` (`Dataset`). 
 68          """ 
 69          # Set up the environment for fitting the data 
 70          X = data.samples.T 
 71          d = data.labels 
 72          if set(d) != set([0, 1]): 
 73              raise ValueError, \ 
 74                    "Regressors for logistic regression should be [0,1]. Got %s" \ 
 75                    %(set(d),) 
 76   
 77          if self.__reduced != 0 : 
 78              # Data have reduced rank 
 79              from scipy.linalg import svd 
 80   
 81              # Compensate for reduced rank: 
 82              # Select only the n largest eigenvectors 
 83              U, S, V = svd(X.T) 
 84              if S[0] == 0: 
 85                  raise FailedToTrainError( 
 86                      "Data provided to PLR seems to be degenerate -- " 
 87                      "0-th singular value is 0") 
 88              S /= S[0] 
 89              V = N.matrix(V[:, :N.max(N.where(S > self.__reduced)) + 1]) 
 90              # Map Data to the subspace spanned by the eigenvectors 
 91              X = (X.T * V).T 
 92   
 93          nfeatures, npatterns = X.shape 
 94   
 95          # Weighting vector 
 96          w  = N.matrix(N.zeros( (nfeatures + 1, 1), 'd')) 
 97          # Error for convergence criterion 
 98          dw = N.matrix(N.ones(  (nfeatures + 1, 1), 'd')) 
 99          # Patterns of interest in the columns 
100          X = N.matrix( \ 
101                  N.concatenate((X, N.ones((1, npatterns), 'd')), 0) \ 
102                  ) 
103          p = N.matrix(N.zeros((1, npatterns), 'd')) 
104          # Matrix implementation of penalty term 
105          Lambda = self.__lm * N.identity(nfeatures + 1, 'd') 
106          Lambda[nfeatures, nfeatures] = 0 
107          # Gradient 
108          g = N.matrix(N.zeros((nfeatures + 1, 1), 'd')) 
109          # Fisher information matrix 
110          H = N.matrix(N.identity(nfeatures + 1, 'd')) 
111   
112          # Optimize 
113          k = 0 
114          while N.sum(N.ravel(dw.A ** 2)) > self.__criterion: 
115              p[:, :] = self.__f(w.T * X) 
116              g[:, :] = X * (d - p).T - Lambda * w 
117              H[:, :] = X * N.diag(p.A1 * (1 - p.A1)) * X.T + Lambda 
118              dw[:, :] = H.I * g 
119              w += dw 
120              k += 1 
121              if k > self.__maxiter: 
122                  raise ConvergenceError, \ 
123                        "More than %d Iterations without convergence" % \ 
124                        (self.__maxiter) 
125   
126          if __debug__: 
127              debug("PLR", \ 
128                    "PLR converged after %d steps. Error: %g" % \ 
129                    (k, N.sum(N.ravel(dw.A ** 2)))) 
130   
131          if self.__reduced: 
132              # We have computed in rank reduced space -> 
133              # Project to original space 
134              self.w = V * w[:-1] 
135              self.offset = w[-1] 
136          else: 
137              self.w = w[:-1] 
138              self.offset = w[-1] 
139   
140   
141 -    def __f(self, y): 
142          """This is the logistic function f, that is used for determination of 
143          the vector w""" 
144          return 1. / (1 + N.exp(-y)) 
145   
146   
147 -    def _predict(self, data): 
148          """ 
149          Predict the class labels for the provided data 
150   
151          Returns a list of class labels 
152          """ 
153          # make sure the data are in matrix form 
154          data = N.matrix(N.asarray(data)) 
155   
156          # get the values and then predictions 
157          values = N.ravel(self.__f(self.offset + data * self.w)) 
158          predictions = values > 0.5 
159   
160          # save the state if desired, relying on State._setitem_ to 
161          # decide if we will actually save the values 
162          self.predictions = predictions 
163          self.values = values 
164   
165          return predictions 
166