Package mvpa :: Package clfs :: Module blr

[hide private]
[frames] | no frames]

# 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