Package mvpa :: Package clfs :: Module gpr
[hide private]
[frames] | no frames]

Source Code for Module mvpa.clfs.gpr

  1  #emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*- 
  2  #ex: set 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  """Gaussian Process Regression (GPR).""" 
 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  from mvpa.clfs.kernel import KernelSquaredExponential 
 20   
 21  if __debug__: 
 22      from mvpa.misc import debug 
 23   
 24   
25 -class GPR(Classifier):
26 """Gaussian Process Regression (GPR). 27 28 """ 29 30 predicted_variances = StateVariable(enabled=False, 31 doc="Variance per each predicted value") 32 33 log_marginal_likelihood = StateVariable(enabled=False, 34 doc="Log Marginal Likelihood") 35 36 37 _clf_internals = [ 'gpr', 'regression', 'non-linear' ] 38
39 - def __init__(self, kernel=KernelSquaredExponential(), 40 sigma_noise=0.001, **kwargs):
41 """Initialize a GPR regression analysis. 42 43 :Parameters: 44 kernel : Kernel 45 a kernel object defining the covariance between instances. 46 (Defaults to KernelSquaredExponential()) 47 sigma_noise : float 48 the standard deviation of the gaussian noise. 49 (Defaults to 0.001) 50 51 """ 52 # init base class first 53 Classifier.__init__(self, **kwargs) 54 55 # pylint happiness 56 self.w = None 57 58 # It does not make sense to calculate a confusion matrix for a GPR 59 self.states.enable('training_confusion', False) 60 61 # set kernel: 62 self.__kernel = kernel 63 64 # set noise level: 65 self.sigma_noise = sigma_noise 66 67 self.predicted_variances = None 68 self.log_marginal_likelihood = None 69 self.train_fv = None 70 self.labels = None 71 self.km_train_train = None 72 pass
73
74 - def __repr__(self):
75 """String summary of the object 76 """ 77 return """GPR(kernel=%s, sigma_noise=%f, enable_states=%s)""" % \ 78 (self.__kernel, self.sigma_noise, str(self.states.enabled))
79 80
82 """ 83 Compute log marginal likelihood using self.train_fv and self.labels. 84 """ 85 # note scipy.cho_factor and scipy.cho_solve seems to be more appropriate 86 # but preliminary tests show them to be slower. 87 88 log_marginal_likelihood = -0.5*N.dot(self.train_labels, self.alpha) - \ 89 N.log(self.L.diagonal()).sum() - \ 90 self.km_train_train.shape[0]*0.5*N.log(2*N.pi) 91 self.log_marginal_likelihood = log_marginal_likelihood 92 93 return log_marginal_likelihood
94 95
96 - def _train(self, data):
97 """Train the classifier using `data` (`Dataset`). 98 """ 99 100 self.train_fv = data.samples 101 self.train_labels = data.labels 102 103 if __debug__: 104 debug("GPR", "Computing train train kernel matrix") 105 106 self.km_train_train = self.__kernel.compute(self.train_fv) 107 108 self.L = N.linalg.cholesky(self.km_train_train + 109 self.sigma_noise**2*N.identity(self.km_train_train.shape[0], 'd')) 110 self.alpha = N.linalg.solve(self.L.transpose(), 111 N.linalg.solve(self.L, self.train_labels)) 112 113 # compute only if the state is enabled 114 if self.states.isEnabled('log_marginal_likelihood'): 115 self.compute_log_marginal_likelihood() 116 117 pass
118 119
120 - def _predict(self, data):
121 """ 122 Predict the output for the provided data. 123 """ 124 if __debug__: 125 debug('GPR', "Computing train test kernel matrix") 126 km_train_test = self.__kernel.compute(self.train_fv, data) 127 if __debug__: 128 debug('GPR', "Computing test test kernel matrix") 129 km_test_test = self.__kernel.compute(data) 130 131 predictions = N.dot(km_train_test.transpose(), self.alpha) 132 133 if self.states.isEnabled('predicted_variances'): 134 # do computation only if state variable was enabled 135 v = N.linalg.solve(self.L, km_train_test) 136 self.predicted_variances = \ 137 N.diag(km_test_test-N.dot(v.transpose(), v)) 138 139 return predictions
140 141
142 - def set_hyperparameters(self,*args):
143 """ 144 Set hyperparameters' values. 145 146 Note that this is a list so the order of the values is 147 important. 148 """ 149 args=args[0] 150 self.sigma_noise = args[0] 151 if len(args)>1: 152 self.__kernel.set_hyperparameters(*args[1:]) 153 pass 154 return
155 156 pass
157 158
159 -def gen_data(n_instances, n_features, flat=False, noise=0.4):
160 """ 161 Generate a (quite) complex multidimensional dataset. 162 """ 163 data = None 164 if flat: 165 data = (N.arange(0.0, 1.0, 1.0/n_instances)*N.pi) 166 data.resize(n_instances, n_features) 167 # print data 168 else: 169 data = N.random.rand(n_instances, n_features)*N.pi 170 pass 171 label = N.sin((data**2).sum(1)).round() 172 label += N.random.rand(label.size)*noise 173 return data, label
174 175
176 -def compute_prediction(sigma_noise_best,length_scale_best,regression,dataset,data_test,label_test,F,logml=True):
177 data_train = dataset.samples 178 label_train = dataset.labels 179 import pylab 180 kse = KernelSquaredExponential(length_scale=length_scale_best) 181 g = GPR(kse, sigma_noise=sigma_noise_best,regression=regression) 182 print g 183 if regression: 184 g.states.enable("predicted_variances") 185 pass 186 187 if logml: 188 g.states.enable("log_marginal_likelihood") 189 pass 190 191 g.train(dataset) 192 prediction = g.predict(data_test) 193 194 # print label_test 195 # print prediction 196 accuracy = None 197 if regression: 198 accuracy = N.sqrt(((prediction-label_test)**2).sum()/prediction.size) 199 print "RMSE:",accuracy 200 else: 201 accuracy = (prediction.astype('l')==label_test.astype('l')).sum()/float(prediction.size) 202 print "accuracy:", accuracy 203 pass 204 205 if F == 1: 206 pylab.figure() 207 pylab.plot(data_train, label_train, "ro", label="train") 208 pylab.plot(data_test, prediction, "b+-", label="prediction") 209 pylab.plot(data_test, label_test, "gx-", label="test (true)") 210 if regression: 211 pylab.plot(data_test, prediction-N.sqrt(g.predicted_variances), "b--", label=None) 212 pylab.plot(data_test, prediction+N.sqrt(g.predicted_variances), "b--", label=None) 213 pylab.text(0.5, -0.8, "RMSE="+"%f" %(accuracy)) 214 else: 215 pylab.text(0.5, -0.8, "accuracy="+str(accuracy)) 216 pass 217 pylab.legend() 218 pass 219 220 print "LML:",g.log_marginal_likelihood
221 222 223 224 225 if __name__ == "__main__": 226 227 N.random.seed(1) 228 229 import pylab 230 pylab.close("all") 231 pylab.ion() 232 233 from mvpa.datasets import Dataset 234 235 train_size = 40 236 test_size = 100 237 F = 1 238 239 data_train, label_train = gen_data(train_size, F) 240 # print label_train 241 242 data_test, label_test = gen_data(test_size, F, flat=True) 243 # print label_test 244 245 dataset = Dataset(samples=data_train, labels=label_train) 246 247 regression = True 248 logml = True 249 250 if logml : 251 print "Looking for better hyperparameters: grid search" 252 253 sigma_noise_steps = N.linspace(0.1, 0.6, num=20) 254 length_scale_steps = N.linspace(0.01, 1.0, num=20) 255 lml = N.zeros((len(sigma_noise_steps), len(length_scale_steps))) 256 lml_best = -N.inf 257 length_scale_best = 0.0 258 sigma_noise_best = 0.0 259 i = 0 260 for x in sigma_noise_steps: 261 j = 0 262 for y in length_scale_steps: 263 kse = KernelSquaredExponential(length_scale=y) 264 g = GPR(kse, sigma_noise=x,regression=regression) 265 g.states.enable("log_marginal_likelihood") 266 g.train(dataset) 267 lml[i,j] = g.log_marginal_likelihood 268 # print x,y,g.log_marginal_likelihood 269 # g.train_fv = dataset.samples 270 # g.train_labels = dataset.labels 271 # lml[i,j] = g.compute_log_marginal_likelihood() 272 if lml[i,j] > lml_best: 273 lml_best = lml[i,j] 274 length_scale_best = y 275 sigma_noise_best = x 276 # print x,y,lml_best 277 pass 278 j += 1 279 pass 280 i += 1 281 pass 282 pylab.figure() 283 X = N.repeat(sigma_noise_steps[:,N.newaxis],sigma_noise_steps.size,axis=1) 284 Y = N.repeat(length_scale_steps[N.newaxis,:],length_scale_steps.size,axis=0) 285 step = (lml.max()-lml.min())/30 286 pylab.contour(X,Y, lml, N.arange(lml.min(), lml.max()+step, step),colors='k') 287 pylab.plot([sigma_noise_best],[length_scale_best],"k+",markeredgewidth=2, markersize=8) 288 pylab.xlabel("noise standard deviation") 289 pylab.ylabel("characteristic length_scale") 290 pylab.title("log marginal likelihood") 291 pylab.axis("tight") 292 print "lml_best",lml_best 293 print "sigma_noise_best",sigma_noise_best 294 print "length_scale_best",length_scale_best 295 print "number of expected upcrossing on the unitary intervale:",1.0/(2*N.pi*length_scale_best) 296 pass 297 298 299 300 compute_prediction(sigma_noise_best,length_scale_best,regression,dataset,data_test,label_test,F,logml) 301