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

Source Code for Module mvpa.clfs.model_selector

  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  """Model selction.""" 
 11   
 12  __docformat__ = 'restructuredtext' 
 13   
 14   
 15  import numpy as N 
 16  from scikits.openopt import NLP 
 17   
 18   
19 -class ModelSelector(object):
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
35 - def max_log_marginal_likelihood(self, hyp_initial_guess,maxiter=1, optimization_algorithm="scipy_cg",ftol=1.0e-2):
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 # minus sign because optimizers do _minimization_.
70 71 def df(*args): 72 """ 73 Proxy to the log_marginal_likelihood first 74 derivative. Necessary for OpenOpt when using derivatives. 75 """ 76 # TODO 77 return
78 79 x0 = hyp_initial_guess # vector of hyperparameters' values where to start the search 80 self.problem = NLP(f,x0) # actual instance of the OpenOpt non-linear problem 81 self.problem.lb = N.zeros(self.problem.n) # set lower bound for hyperparameters: avoid negative hyperparameters. Note: problem.n is the size of hyperparameters' vector 82 self.problem.maxiter = maxiter # max number of iterations for the optimizer. 83 self.problem.checkdf = True # check whether the derivative of log_marginal_likelihood converged to zero before ending optimization 84 self.problem.ftol = ftol # set increment of log_marginal_likelihood under which the optimizer stops 85 self.problem.iprint = 0 # shut up OpenOpt (note: -1 = no logs, 0 = small log, 1 = verbose) 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) # perform optimization! 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 # best hyperparameters found # NOTE is it better to return a copy? 99 self.log_marginal_likelihood_best = -result.ff # actual best vuale of log_marginal_likelihood 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 # print label_train 128 129 data_test, label_test = gpr.gen_data(test_size, F, flat=True) 130 # print label_test 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 # k = kernel.KernelConstant() 139 g = gpr.GPR(k,regression=regression) 140 g.states.enable("log_marginal_likelihood") 141 # g.train_fv = dataset.samples 142 # g.train_labels = dataset.labels 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 # data = N.hstack([data]*10) # test a larger set of dimensions: reduce ftol! 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