Source Code for Module mvpa.clfs.blr

  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  #   Copyright (c) 2008 Emanuele Olivetti <emanuele@relativita.com> 
  6  #   See COPYING file distributed along with the PyMVPA package for the 
  7  #   copyright and license terms. 
  8  # 
  9  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 
 10  """Bayesian Linear Regression (BLR).""" 
 11   
 12  __docformat__ = 'restructuredtext' 
 13   
 14   
 15  import numpy as N 
 16   
 17  from mvpa.misc.state import StateVariable 
 18  from mvpa.clfs.base import Classifier 
 19   
 20  if __debug__: 
 21      from mvpa.misc import debug 
 22   
 23   
 24 -class BLR(Classifier): 
 25      """Bayesian Linear Regression (BLR). 
 26   
 27      """ 
 28   
 29      predicted_variances = StateVariable(enabled=False, 
 30          doc="Variance per each predicted value") 
 31   
 32      log_marginal_likelihood = StateVariable(enabled=False, 
 33          doc="Log Marginal Likelihood") 
 34   
 35   
 36      _clf_internals = [ 'blr', 'regression', 'linear' ] 
 37   
 38 -    def __init__(self, sigma_p = None, sigma_noise=1.0, **kwargs): 
 39          """Initialize a BLR regression analysis. 
 40   
 41          :Parameters: 
 42            sigma_noise : float 
 43              the standard deviation of the gaussian noise. 
 44              (Defaults to 0.1) 
 45   
 46          """ 
 47          # init base class first 
 48          Classifier.__init__(self, **kwargs) 
 49   
 50          # pylint happiness 
 51          self.w = None 
 52   
 53          # It does not make sense to calculate a confusion matrix for a 
 54          # BLR: 
 55          self.states.enable('training_confusion', False) 
 56   
 57          # set the prior on w: N(0,sigma_p) , specifying the covariance 
 58          # sigma_p on w: 
 59          self.sigma_p = sigma_p 
 60   
 61          # set noise level: 
 62          self.sigma_noise = sigma_noise 
 63   
 64          self.predicted_variances = None 
 65          self.log_marginal_likelihood = None 
 66          self.labels = None 
 67          pass 
 68   
 69 -    def __repr__(self): 
 70          """String summary of the object 
 71          """ 
 72          return """BLR(w=%s, sigma_p=%s, sigma_noise=%f, enable_states=%s)""" % \ 
 73                 (self.w, self.sigma_p, self.sigma_noise, str(self.states.enabled)) 
 74   
 75   
 76 -    def compute_log_marginal_likelihood(self): 
 77          """ 
 78          Compute log marginal likelihood using self.train_fv and self.labels. 
 79          """ 
 80          # log_marginal_likelihood = None 
 81          # return log_marginal_likelihood 
 82          raise NotImplementedError 
 83       
 84   
 85 -    def _train(self, data): 
 86          """Train regression using `data` (`Dataset`). 
 87          """ 
 88          # provide a basic (i.e. identity matrix) and correct prior 
 89          # sigma_p, if not provided before or not compliant to 'data': 
 90          if self.sigma_p == None: # case: not provided 
 91              self.sigma_p = N.eye(data.samples.shape[1]+1) 
 92          elif self.sigma_p.shape[1] != (data.samples.shape[1]+1): # case: wrong dimensions 
 93              self.sigma_p = N.eye(data.samples.shape[1]+1) 
 94          else: 
 95              # ...then everything is OK :) 
 96              pass 
 97   
 98          # add one fake column of '1.0' to model the intercept: 
 99          self.samples_train = N.hstack([data.samples,N.ones((data.samples.shape[0],1))]) 
100          if type(self.sigma_p)!=type(self.samples_train): # if sigma_p is a number... 
101              self.sigma_p = N.eye(self.samples_train.shape[1])*self.sigma_p # convert in matrix 
102              pass 
103   
104          self.A_inv = N.linalg.inv(1.0/(self.sigma_noise**2) * 
105                                    N.dot(self.samples_train.T, 
106                                          self.samples_train) + 
107                                    N.linalg.inv(self.sigma_p)) 
108          self.w = 1.0/(self.sigma_noise**2) * N.dot(self.A_inv, 
109                                                     N.dot(self.samples_train.T, 
110                                                           data.labels)) 
111          pass 
112   
113   
114 -    def _predict(self, data): 
115          """ 
116          Predict the output for the provided data. 
117          """ 
118   
119          data = N.hstack([data,N.ones((data.shape[0],1),dtype=data.dtype)]) 
120          predictions = N.dot(data,self.w) 
121           
122          if self.states.isEnabled('predicted_variances'): 
123              # do computation only if state variable was enabled 
124              self.predicted_variances = N.dot(data, N.dot(self.A_inv, data.T)).diagonal()[:,N.newaxis] 
125   
126          return predictions 
127   
128   
129 -    def set_hyperparameters(self,*args): 
130          """ 
131          Set hyperparameters' values. 
132   
133          Note that this is a list so the order of the values is 
134          important. 
135          """ 
136          args=args[0] 
137          self.sigma_noise = args[0] 
138          if len(args)>1: 
139              self.sigma_p = N.array(args[1:]) # XXX check if this is ok 
140              pass 
141          return 
142   
143      pass 
144   
145   
146  if __name__ == "__main__": 
147      import pylab 
148      pylab.close("all") 
149      pylab.ion() 
150   
151      from mvpa.misc.data_generators import linear_awgn 
152   
153      train_size = 10 
154      test_size = 100 
155      F = 1 # dimensions of the dataset 
156   
157      # N.random.seed(1) 
158   
159      slope = N.random.rand(F) 
160      intercept = N.random.rand(1) 
161      print "True slope:",slope 
162      print "True intercept:",intercept 
163       
164      dataset_train = linear_awgn(train_size, intercept=intercept, slope=slope) 
165      # print dataset.labels 
166   
167      dataset_test = linear_awgn(test_size, intercept=intercept, slope=slope, flat=True) 
168   
169      regression = True 
170      logml = False 
171   
172      b = BLR(sigma_p=N.eye(F+1), sigma_noise=0.1, regression=True) 
173      b.states.enable("predicted_variances") 
174      b.train(dataset_train) 
175      predictions = b.predict(dataset_test.samples) 
176      print "Predicted slope and intercept:",b.w 
177       
178      if F==1: 
179          pylab.plot(dataset_train.samples,dataset_train.labels,"ro",label="train") 
180           
181          pylab.plot(dataset_test.samples,predictions,"b-",label="prediction") 
182          pylab.plot(dataset_test.samples,predictions+N.sqrt(b.predicted_variances),"b--",label="pred(+/-)std") 
183          pylab.plot(dataset_test.samples,predictions-N.sqrt(b.predicted_variances),"b--",label=None) 
184          # pylab.plot(dataset_test.samples,dataset_test.labels,"go") 
185          pylab.legend() 
186          pylab.xlabel("samples") 
187          pylab.ylabel("labels") 
188          pylab.title("Bayesian Linear Regression on dataset 'linear_AWGN'") 
189          pass 
190