1
2
3
4
5
6
7
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
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
53 Classifier.__init__(self, **kwargs)
54
55
56 self.w = None
57
58
59 self.states.enable('training_confusion', False)
60
61
62 self.__kernel = kernel
63
64
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
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
94
95
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
114 if self.states.isEnabled('log_marginal_likelihood'):
115 self.compute_log_marginal_likelihood()
116
117 pass
118
119
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
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
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
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
195
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
241
242 data_test, label_test = gen_data(test_size, F, flat=True)
243
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
269
270
271
272 if lml[i,j] > lml_best:
273 lml_best = lml[i,j]
274 length_scale_best = y
275 sigma_noise_best = x
276
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