


Run Gaussian process model - wrapper for gpml toolbox
FORMAT output = prt_machine_gpml(d,args)
Inputs:
d - structure with data information, with mandatory fields:
.train - training data (cell array of matrices of row vectors,
each [Ntr x D]). each matrix contains one representation
of the data. This is useful for approaches such as
multiple kernel learning.
.test - testing data (cell array of matrices row vectors, each
[Nte x D])
.testcov - testing covariance (cell array of matrices row vectors,
each [Nte x Nte])
.tr_targets - training labels (for classification) or values (for
regression) (column vector, [Ntr x 1])
.use_kernel - flag, is data in form of kernel matrices (true) or in
form of features (false)
args - argument string, where
-h - optimise hyperparameters (otherwise don't)
-f iter - max # iterations for optimiser (ignored if -h not set)
-l likfun - likelihood function:
'likErf' - erf/probit likelihood (binary only)
-c covfun - covariance function:
'covLINkcell' - simple dot product
'covLINglm' - construct a GLM
-m meanfun - mean function:
'meanConstcell' - suitable for dot product
'meanConstglm' - suitable for GLM
-i inffun - inference function:
'prt_infEP' - Expectation Propagation
experimental args (use at your own risk):
-p - use priors for the hyperparameters. If specified, this
indicates that a maximum a posteriori (MAP) approach
will be used to set covariance function
hyperparameters. The priors are obtained by calling
prt_gp_priors('covFuncName')
N.B.: for the arguments specifying functions, pass in a string, not
a function handle. This script will generate a function handle
Output:
output - output of machine (struct).
* Mandatory fields:
.predictions - predictions of classification or regression [Nte x D]
* Optional fields:
.type - which type of machine this is (here, 'classifier')
.func_val - predictive probabilties
.mu - test latent means
.s2 - test latent variances
.loghyper - log hyperparameters
.nlml - negative log marginal likelihood
.alpha - GP weighting coefficients
.sW - likelihood matrix (see Rasmussen & Williams, 2006)
.L - Cholesky factor
__________________________________________________________________________
Copyright (C) 2011 Machine Learning & Neuroimaging Laboratory

0001 function output = prt_machine_gpml(d,args) 0002 % Run Gaussian process model - wrapper for gpml toolbox 0003 % FORMAT output = prt_machine_gpml(d,args) 0004 % Inputs: 0005 % d - structure with data information, with mandatory fields: 0006 % .train - training data (cell array of matrices of row vectors, 0007 % each [Ntr x D]). each matrix contains one representation 0008 % of the data. This is useful for approaches such as 0009 % multiple kernel learning. 0010 % .test - testing data (cell array of matrices row vectors, each 0011 % [Nte x D]) 0012 % .testcov - testing covariance (cell array of matrices row vectors, 0013 % each [Nte x Nte]) 0014 % .tr_targets - training labels (for classification) or values (for 0015 % regression) (column vector, [Ntr x 1]) 0016 % .use_kernel - flag, is data in form of kernel matrices (true) or in 0017 % form of features (false) 0018 % args - argument string, where 0019 % -h - optimise hyperparameters (otherwise don't) 0020 % -f iter - max # iterations for optimiser (ignored if -h not set) 0021 % -l likfun - likelihood function: 0022 % 'likErf' - erf/probit likelihood (binary only) 0023 % -c covfun - covariance function: 0024 % 'covLINkcell' - simple dot product 0025 % 'covLINglm' - construct a GLM 0026 % -m meanfun - mean function: 0027 % 'meanConstcell' - suitable for dot product 0028 % 'meanConstglm' - suitable for GLM 0029 % -i inffun - inference function: 0030 % 'prt_infEP' - Expectation Propagation 0031 % experimental args (use at your own risk): 0032 % -p - use priors for the hyperparameters. If specified, this 0033 % indicates that a maximum a posteriori (MAP) approach 0034 % will be used to set covariance function 0035 % hyperparameters. The priors are obtained by calling 0036 % prt_gp_priors('covFuncName') 0037 % 0038 % N.B.: for the arguments specifying functions, pass in a string, not 0039 % a function handle. This script will generate a function handle 0040 % 0041 % Output: 0042 % output - output of machine (struct). 0043 % * Mandatory fields: 0044 % .predictions - predictions of classification or regression [Nte x D] 0045 % * Optional fields: 0046 % .type - which type of machine this is (here, 'classifier') 0047 % .func_val - predictive probabilties 0048 % .mu - test latent means 0049 % .s2 - test latent variances 0050 % .loghyper - log hyperparameters 0051 % .nlml - negative log marginal likelihood 0052 % .alpha - GP weighting coefficients 0053 % .sW - likelihood matrix (see Rasmussen & Williams, 2006) 0054 % .L - Cholesky factor 0055 %__________________________________________________________________________ 0056 % Copyright (C) 2011 Machine Learning & Neuroimaging Laboratory 0057 0058 % Written by A Marquand 0059 % $Id$ 0060 0061 % configure default parameters for GP optimisation 0062 meanfunc = @meanConstcell; 0063 covfunc = @covLINkcell; 0064 maxeval = -20; 0065 if strcmp(d.pred_type,'classification') 0066 mode = 'classifier'; % it's good to be consistent! 0067 likfunc = @likErf; 0068 inffunc = @prt_infEP; 0069 else 0070 mode = 'regression'; 0071 likfunc = @likGauss; 0072 inffunc = @prt_infExact; 0073 mtr = mean(d.tr_targets); % mean of the training data 0074 end 0075 0076 % Error checks 0077 % ------------------------------------------------------------------------- 0078 SANITYCHECK=true; % can turn off for "speed". Expert only. 0079 0080 if SANITYCHECK==true 0081 % args should be a string (empty or otherwise) 0082 if ~ischar(args) 0083 error('prt_machine_gpml:argsNotString',['Error: gpml'... 0084 ' args should be a string. ' ... 0085 ' SOLUTION: Please do XXX']); 0086 end 0087 0088 % check we can reach the binary library 0089 if ~exist('prt_gp','file') 0090 error('prt_machine_gpml:libNotFound',['Error:'... 0091 ' ''prt_gp'' function could not be found !' ... 0092 ' SOLUTION: Please check your path.']); 0093 end 0094 % check whether it is a two-class classification problem 0095 uTL=unique(d.tr_targets(:)); 0096 k=numel(uTL); % number of classes 0097 if strcmp(mode,'classifier') && k > 2 0098 warning('prt_machine_gpml:classificationWithMoreThanTwoClasses',... 0099 ['Classification specified with > 2 classes. ',... 0100 'Defaulting to multiclass Laplace approximation.']); 0101 output = prt_machine_gpclap(d,args); 0102 return; 0103 end 0104 % are we using a kernel ? 0105 if ~d.use_kernel 0106 error('prt_machine_gpml:useKernelIsFalse',['Error:'... 0107 ' This machine is currently only implemented for kernel data ' ... 0108 'SOLUTION: Please set use_kernel to true']); 0109 end 0110 end 0111 0112 % parse input arguments 0113 % ------------------------------------------------------------------------- 0114 % hyperparameters 0115 if ~isempty(regexp(args,'-h','once')) 0116 optimise_theta = true; 0117 eargs = regexp(args,'-f\s+[0-9]*','match'); 0118 if ~isempty(eargs) 0119 eargs = regexp(cell2mat(eargs),'-f\s+','split'); 0120 maxeval = str2num(['-',cell2mat(eargs(2))]); 0121 end 0122 else 0123 optimise_theta = false; 0124 end 0125 % likelihood function 0126 largs = regexp(args,'-l\s+[a-zA-Z0-9_]*','match'); 0127 if ~isempty(largs) 0128 largs = regexp(cell2mat(largs),'-l\s+','split'); 0129 likfunc = str2func(cell2mat(largs(2))); 0130 if strcmpi(cell2mat(largs(2)),'Erf') 0131 likfunc = @likErf; 0132 end 0133 end 0134 % covariance function 0135 cargs = regexp(args,'-c\s+[a-zA-Z0-9_]*','match'); 0136 if ~isempty(cargs) 0137 cargs = regexp(cell2mat(cargs),'-c\s+','split'); 0138 covfunc = str2func(cell2mat(cargs(2))); 0139 end 0140 % mean function 0141 margs = regexp(args,'-m\s+[a-zA-Z0-9_]*','match'); 0142 if ~isempty(margs) 0143 margs = regexp(cell2mat(margs),'-m\s+','split'); 0144 meanfunc = str2func(cell2mat(margs(2))); 0145 end 0146 % inference function 0147 iargs = regexp(args,'-i\s+[a-zA-Z0-9_]*','match'); 0148 if ~isempty(iargs) 0149 iargs = regexp(cell2mat(iargs),'-i\s+','split'); 0150 inffunc = str2func(cell2mat(iargs(2))); 0151 end 0152 % priors 0153 if ~isempty(regexp(args,'-p','once')) 0154 disp('Empirical priors specified. Using MAP for hyperparameters') 0155 priors = prt_gp_priors(func2str(covfunc)); 0156 map = true; 0157 else 0158 map = false; 0159 end 0160 0161 % Set default hyperparameters 0162 % ------------------------------------------------------------------------- 0163 nhyp = str2num([feval(covfunc); feval(likfunc); feval(meanfunc)]); 0164 if nhyp(1) > 0 0165 hyp.cov = zeros(nhyp(1),1); 0166 end 0167 if nhyp(2) > 0 0168 hyp.lik = zeros(nhyp(2),1); 0169 end 0170 if nhyp(3) > 0 0171 hyp.mean = zeros(nhyp(3),1); 0172 end 0173 0174 % Assemble data matrices 0175 % ------------------------------------------------------------------------- 0176 % handle the glm as a special case (for now) 0177 if strcmpi(func2str(covfunc),'covLINglm') || strcmpi(func2str(covfunc),'covLINglm_2class') 0178 % configure covariances 0179 K = [d.train(:)' {d.tr_param}]; 0180 Ks = [d.test(:)' {d.te_param}]; 0181 Kss = [d.testcov(:)' {d.te_param}]; 0182 0183 % get default hyperparamter values 0184 hyp.cov = log(prt_glm_design); 0185 0186 [tmp1 tmp2 tmp3 tr_lbs] = prt_glm_design(hyp.cov, d.tr_param); 0187 [tmp1 tmp2 tmp3 te_lbs] = prt_glm_design(hyp.cov, d.te_param); 0188 else 0189 % configure covariances 0190 K = d.train; 0191 Ks = d.test; 0192 Kss = d.testcov; 0193 0194 tr_lbs = d.tr_targets; 0195 te_lbs = d.te_targets; 0196 end 0197 0198 % configure targets 0199 if strcmp(mode,'classifier') 0200 % convert targets to +1/-1 0201 y = -1*(2 * tr_lbs - 3); 0202 else 0203 y = tr_lbs - mtr; 0204 end 0205 0206 % Train and test GP model 0207 % ------------------------------------------------------------------------- 0208 % train 0209 if optimise_theta 0210 if map 0211 [hyp,nlmls] = minimize(hyp, @prt_gp_map, maxeval, inffunc, meanfunc, covfunc, likfunc, K, y, priors); 0212 else 0213 [hyp nlmls] = minimize(hyp, @prt_gp, maxeval, inffunc, meanfunc, covfunc, likfunc, K, y); 0214 end 0215 else 0216 nlmls = prt_gp(hyp, inffunc, meanfunc, covfunc, likfunc, K, y); 0217 end 0218 0219 % make predictions 0220 [ymu ys2 fmu fs2 lp post] = prt_gp(hyp, inffunc, meanfunc, covfunc, likfunc,K, y, Ks, zeros(size(Ks{1},1),1), Kss); 0221 0222 % Outputs 0223 % ------------------------------------------------------------------------- 0224 if strcmp(mode,'classifier') 0225 p = exp(lp); 0226 output.predictions = (1-real(p > 0.5)) + 1; 0227 output.func_val = p; 0228 else % regression 0229 output.predictions = ymu + mtr; 0230 output.func_val = output.predictions; 0231 end 0232 output.type = mode; 0233 output.loghyper = hyp; 0234 output.mu = ymu; 0235 output.s2 = ys2; 0236 output.nlml = min(nlmls); 0237 output.tr_targets = tr_lbs; 0238 output.te_targets = te_lbs; 0239 output.alpha = post.alpha; 0240 %output.sW = post.sW; 0241 %output.L = post.L; 0242 0243 end 0244