function result=cosmo_naive_bayes_classifier_searchlight(ds, nbrhood, varargin)
% Run (fast) Naive Bayes classifier searchlight with crossvalidation
%
% result=cosmo_naive_bayes_classifier_searchlight(ds, nbrhood, ...)
%
% Inputs:
% ds dataset struct
% nbrhood Neighborhood structure with fields:
% .a struct with dataset attributes
% .fa struct with feature attributes. Each field
% should have NF values in the second dimension
% .neighbors cell with NF mappings from center_ids in output
% dataset to feature ids in input dataset.
% Suitable neighborhood structs can be generated
% using:
% - cosmo_spherical_neighborhood (fmri volume)
% - cosmo_surficial_neighborhood (fmri surface)
% - cosmo_meeg_chan_neigborhood (MEEG channels)
% - cosmo_interval_neighborhood (MEEG time, freq)
% - cosmo_cross_neighborhood (to cross neighborhoods
% from the neighborhood
% functions above)
% 'partitions', par Partition scheme to use. Typically this is the
% output from cosmo_nfold_partitioner or
% cosmo_oddeven_partitioner. Partitions schemes
% with more than one prediction for samples in the
% test set (such as the output from
% cosmo_nchoosek_partitioner(N) with N>1) are not
% supported
% 'output', out One of:
% 'accuracy' return classification accuracy
% 'predictions' return prediction for each sample
% in a test set in partitions
% 'progress', p Show progress every p folds (default: 1)
%
% Output:
% results_map a dataset struct where the samples
% contain classification accuracies or class
% predictions
%
% Example:
% % generate tiny dataset (6 voxels) and define a tiny spherical
% % neighborhood with a radius of 1 voxel.
% ds=cosmo_synthetic_dataset('nchunks',10,'ntargets',5);
% nh=cosmo_spherical_neighborhood(ds,'radius',1,'progress',false);
% %
% % set options
% opt=struct();
% opt.progress=false;
% % define take-one-chunk-out crossvalidation scheme (10 folds)
% opt.partitions=cosmo_nfold_partitioner(ds);
% %
% % run searchlight
% result=cosmo_naive_bayes_classifier_searchlight(ds,nh,opt);
% %
% % show result
% cosmo_disp(result);
% %|| .a
% %|| .fdim
% %|| .labels
% %|| { 'i' 'j' 'k' }
% %|| .values
% %|| { [ 1 2 3 ] [ 1 2 ] [ 1 ] }
% %|| .vol
% %|| .mat
% %|| [ 2 0 0 -3
% %|| 0 2 0 -3
% %|| 0 0 2 -3
% %|| 0 0 0 1 ]
% %|| .dim
% %|| [ 3 2 1 ]
% %|| .xform
% %|| 'scanner_anat'
% %|| .fa
% %|| .nvoxels
% %|| [ 3 4 3 3 4 3 ]
% %|| .radius
% %|| [ 1 1 1 1 1 1 ]
% %|| .center_ids
% %|| [ 1 2 3 4 5 6 ]
% %|| .i
% %|| [ 1 2 3 1 2 3 ]
% %|| .j
% %|| [ 1 1 1 2 2 2 ]
% %|| .k
% %|| [ 1 1 1 1 1 1 ]
% %|| .samples
% %|| [ 0.6 0.74 0.44 0.6 0.7 0.4 ]
% %|| .sa
% %|| .labels
% %|| { 'accuracy' }
%
%
% Notes:
% - this function runs considerably faster than using a searchlight with
% a classifier function and a crossvalidation scheme, because model
% parameters during training are estimated only once for each feature.
% Thus, speedups are most significant if elements in the neighborhood
% have many overlapping features
% - for other classifiers or other measures, use the more flexible
% cosmo_searchlight function
%
% See also: cosmo_searchlight
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
% check input
cosmo_check_dataset(ds);
% set defaults
defaults=struct();
defaults.output='accuracy';
defaults.progress=1;
opt=cosmo_structjoin(defaults,varargin{:});
show_progress=isfield(opt,'progress') && opt.progress;
% get neighborhood in matrix form for faster lookups
nbrhood_mat=cosmo_convert_neighborhood(nbrhood,'matrix');
ncenters=size(nbrhood_mat,2);
nsamples=size(ds.samples,1);
% get partitions for crossvalidation
[train_idxs,test_idxs]=get_partitions(ds,opt);
nfolds=numel(train_idxs);
max_prediction_count=get_max_prediction_count(test_idxs);
predictions=NaN(nsamples,ncenters,max_prediction_count);
prediction_count=zeros(nsamples,1);
if show_progress
clock_start=clock();
prev_progress_msg='';
end
% perform classification for each fold
for fold=1:nfolds
train_idx=train_idxs{fold};
test_idx=test_idxs{fold};
samples_train=ds.samples(train_idx,:);
targets_train=ds.sa.targets(train_idx);
% estimate parameters
model=naive_bayes_train(samples_train, targets_train);
% predict classes
test_samples=ds.samples(test_idx,:);
fold_pred=naive_bayes_predict(model, nbrhood_mat, test_samples);
% store predictions; work backwards to ensure each is stored just
% once
for col=max_prediction_count:-1:1
row_msk=prediction_count(test_idx)==(col-1);
row=test_idx(row_msk);
predictions(row,:,col)=fold_pred(row_msk,:);
prediction_count(row)=prediction_count(row)+1;
end
if show_progress && (fold<10 || ~mod(fold,opt.progress) || ...
fold==ncenters)
msg='';
prev_progress_msg=cosmo_show_progress(clock_start, ...
fold/nfolds, msg, prev_progress_msg);
end
end
result=struct();
result.a=nbrhood.a;
result.fa=nbrhood.fa;
% set output
output=opt.output;
switch output
case 'accuracy'
is_pred=~isnan(predictions);
is_correct=is_pred & bsxfun(@eq,predictions,ds.sa.targets);
correct_count=sum(sum(is_correct,1),3);
pred_count=sum(sum(is_pred,1),3);
result.samples=correct_count./pred_count;
result.sa.labels={'accuracy'};
case {'winner_predictions','predictions'}
if cosmo_match({output},{'predictions'})
cosmo_warning('CoSMoMVPA:deprecated',...
sprintf(...
['Output option ''%s'' is deprecated and will '...
'be removed from a future release. Please use '...
'output=''winner_predictions'' instead.'],...
output));
end
if max_prediction_count<=1
winners=predictions;
else
winners=zeros(nsamples,ncenters);
for k=1:ncenters
pred_mat=reshape(predictions(:,k,:),...
nsamples,max_prediction_count);
[idx,cl]=cosmo_winner_indices(pred_mat);
winners(:,k)=cl(idx);
if show_progress && (k<10 || ...
~mod(k,opt.progress) || ...
k==ncenters)
msg='computing winners';
prev_progress_msg=cosmo_show_progress(...
clock_start, ...
k/ncenters, msg, prev_progress_msg);
end
end
end
result.samples=winners;
result.sa=rmfield(ds.sa,'chunks');
otherwise
error(['illegal output ''%s'', must be '...
'''accuracy'' or ''winner_predictions'''], opt.output);
end
cosmo_check_dataset(result);
function [train_idxs,test_idxs]=get_partitions(ds,opt)
if ~isfield(opt,'partitions')
error('the ''partitions'' option is required');
end
partitions=opt.partitions;
cosmo_check_partitions(partitions,ds);
train_idxs=partitions.train_indices;
test_idxs=partitions.test_indices;
function max_prediction_count=get_max_prediction_count(test_idxs)
max_index=max(cellfun(@max,test_idxs));
h=zeros(max_index,1);
nfolds=numel(test_idxs);
for fold=1:nfolds
idx=test_idxs{fold};
h(idx)=h(idx)+1;
end
max_prediction_count=max(h);
function pred=naive_bayes_predict(model,nbrhood_mat,samples_train)
classes=model.classes;
nclasses=numel(classes);
[max_neighbors,ncenters]=size(nbrhood_mat);
[nsamples,nfeatures]=size(samples_train);
if nfeatures~=size(model.mus,2)
error(['size mismatch in number of features between training '...
'and test set']);
end
max_ps=-Inf(nsamples,ncenters);
pred=NaN(nsamples,ncenters);
for j=1:nclasses
mu=model.mus(j,:);
var_=model.vars(j,:);
% nsamples x nfeatures
xs_z=bsxfun(@rdivide,bsxfun(@minus,samples_train,mu).^2,var_);
log_ps=-.5*(bsxfun(@plus,log(2*pi*var_),xs_z)) + ...
model.log_class_probs(j);
log_sum_ps=zeros(nsamples,ncenters);
for k=1:max_neighbors
row_msk=nbrhood_mat(k,:)>0;
log_sum_ps(:,row_msk)=log_sum_ps(:,row_msk)+...
log_ps(:,nbrhood_mat(k,row_msk));
end
greater_ps_mask=log_sum_ps>max_ps;
max_ps(greater_ps_mask)=log_sum_ps(greater_ps_mask);
pred(greater_ps_mask)=classes(j);
end
function model=naive_bayes_train(samples_train, targets_train)
[ntrain,nfeatures]=size(samples_train);
if ntrain~=numel(targets_train)
error('size mismatch between samples and targets');
end
[class_idxs,classes_cell]=cosmo_index_unique({targets_train});
classes=classes_cell{1};
nclasses=numel(classes);
% allocate space for statistics of each class
mus=zeros(nclasses,nfeatures);
vars=zeros(nclasses,nfeatures);
log_class_probs=zeros(nclasses,1);
% compute means and standard deviations of each class
for k=1:nclasses
idx=class_idxs{k};
nsamples_in_class=numel(idx); % number of samples
if nsamples_in_class<2
error(['Cannot train: class %d has only %d samples, %d '...
'are required'],nsamples_in_class,classes(k));
end
d=samples_train(idx,:); % samples in this class
mu=mean(d); %mean
mus(k,:)=mu;
% variance - faster than 'var'
vars(k,:)=sum(bsxfun(@minus,mu,d).^2,1)/nsamples_in_class;
% log of class probability
log_class_probs(k)=log(nsamples_in_class/ntrain);
end
model.mus=mus;
model.vars=vars;
model.log_class_probs=log_class_probs;
model.classes=classes;
