1
2
3
4
5
6
7
8
9
10 """Model selction."""
11
12 __docformat__ = 'restructuredtext'
13
14
15 import numpy as N
16 from scikits.openopt import NLP
17
18
20 """Model selection facility.
21
22 Select a model among multiple models (i.e., a parametric model,
23 parametrized by a set of hyperparamenters).
24 """
25
26 - def __init__(self, parametric_model, dataset):
27 self.parametric_model = parametric_model
28 self.dataset = dataset
29 self.hyperparameters_best = None
30 self.log_marginal_likelihood_best = None
31 self.problem = None
32 pass
33
34
36 """
37 Set up the optimization problem in order to maximize
38 the log_marginal_likelihood.
39
40 :Parameters:
41
42 parametric_model : Classifier
43 the actual parameteric model to be optimized.
44
45 hyp_initial_guess : numpy.ndarray
46 set of hyperparameters' initial values where to start optimization.
47
48 optimization_algorithm : string
49 actual name of the optimization algorithm. See
50 http://scipy.org/scipy/scikits/wiki/NLP
51 for a comprehensive/updated list of available NLP solvers.
52 (Defaults to 'ralg')
53
54 NOTE: the maximization of log_marginal_likelihood is a non-linear
55 optimization problem (NLP). This fact is confirmed by Dmitrey,
56 author of OpenOpt.
57 """
58 self.optimization_algorithm = optimization_algorithm
59
60 def f(*args):
61 """
62 Wrapper to the log_marginal_likelihood to be
63 maximized. Necessary for OpenOpt since it performs
64 minimization only.
65 """
66 self.parametric_model.set_hyperparameters(*args)
67 self.parametric_model.train(self.dataset)
68 log_marginal_likelihood = self.parametric_model.compute_log_marginal_likelihood()
69 return -log_marginal_likelihood
70
71 def df(*args):
72 """
73 Proxy to the log_marginal_likelihood first
74 derivative. Necessary for OpenOpt when using derivatives.
75 """
76
77 return
78
79 x0 = hyp_initial_guess
80 self.problem = NLP(f,x0)
81 self.problem.lb = N.zeros(self.problem.n)
82 self.problem.maxiter = maxiter
83 self.problem.checkdf = True
84 self.problem.ftol = ftol
85 self.problem.iprint = 0
86 return self.problem
87
88
89 - def solve(self, problem=None):
90 """Solve the minimization problem, check outcome and collect results.
91 """
92 result = self.problem.solve(self.optimization_algorithm)
93 if result.stopcase == -1:
94 print "Unable to find a maximum to log_marginal_likelihood"
95 elif result.stopcase == 0:
96 print "Limits exceeded"
97 elif result.stopcase == 1:
98 self.hyperparameters_best = result.xf
99 self.log_marginal_likelihood_best = -result.ff
100
101 return self.log_marginal_likelihood_best
102
103 pass
104
105
106
107 if __name__ == "__main__":
108
109 import gpr
110 import kernel
111
112 N.random.seed(1)
113
114 import pylab
115 pylab.close("all")
116 pylab.ion()
117
118 from mvpa.datasets import Dataset
119
120 print "GPR:",
121
122 train_size = 40
123 test_size = 100
124 F = 1
125
126 data_train, label_train = gpr.gen_data(train_size, F)
127
128
129 data_test, label_test = gpr.gen_data(test_size, F, flat=True)
130
131
132 dataset = Dataset(samples=data_train, labels=label_train)
133
134 regression = True
135 logml = True
136
137 k = kernel.KernelLinear(coefficient=N.ones(1))
138
139 g = gpr.GPR(k,regression=regression)
140 g.states.enable("log_marginal_likelihood")
141
142
143
144 print "GPR hyperparameters' search through maximization of marginal likelihood on train data."
145 print
146 ms = ModelSelector(g,dataset)
147
148 sigma_noise_initial = 1.0
149 length_scale_initial = 1.0
150
151 problem = ms.max_log_marginal_likelihood(hyp_initial_guess=[sigma_noise_initial,length_scale_initial], optimization_algorithm="ralg", ftol=1.0e-4)
152 lml = ms.solve()
153 sigma_noise_best, length_scale_best = ms.hyperparameters_best
154 print
155 print "Best sigma_noise:",sigma_noise_best
156 print "Best length_scale:",length_scale_best
157 print "Best log_marginal_likelihood:",lml
158
159 gpr.compute_prediction(sigma_noise_best,length_scale_best,regression,dataset,data_test,label_test,F)
160
161 print
162 print "GPR ARD on dataset from Williams and Rasmussen 1996:"
163 data, labels = kernel.generate_dataset_wr1996()
164
165 dataset = Dataset(samples=data, labels=labels)
166 k = kernel.KernelSquaredExponential(length_scale=N.ones(dataset.samples.shape[1]))
167 g = gpr.GPR(k, regression=regression)
168 ms = ModelSelector(g, dataset)
169
170 sigma_noise_initial = 0.01
171 length_scales_initial = 0.5*N.ones(dataset.samples.shape[1])
172
173 problem = ms.max_log_marginal_likelihood(hyp_initial_guess=N.hstack([sigma_noise_initial,length_scales_initial]), optimization_algorithm="ralg")
174 lml = ms.solve()
175 sigma_noise_best = ms.hyperparameters_best[0]
176 length_scales_best = ms.hyperparameters_best[1:]
177 print
178 print "Best sigma_noise:",sigma_noise_best
179 print "Best length_scale:",length_scales_best
180 print "Best log_marginal_likelihood:",lml
181