aic_mdl

(Toolbox/aic_mdl.m in BrainStorm 2.0 (Alpha))

Function Synopsis

[dim, aic, mdl] = aic_mdl(s,n);

Help Text

AIC_MDL - Calculate Information Theoretic Criteria for data dimensionality
 function [dim, aic, mdl] = aic_mdl(s,n);
 Let s be the ordered eigenvalues of the data correlation matrix, calculated from
  n samples of data.
 E.G.: F  is 151 x 375, then s = svd(F).^2; n = 375;
 Alternatively, call as [dim, aic, mdl] = aic_mdl(F);
  where F is the spatio-temporal data matrix.
 Returns the subspace dimension, calculated by two different methods,
  one based on Akaike (AIC), the other by Schwartz and Rissanen (MDL).
 dim is 1 x 2, with the subspace dimension as found by each method
 aic and mdl are the calculated metrics

 See IEEE ASSP 1985, Wax and Kailath.

Cross-Reference Information

parametric_estimation              C:\BrainStorm_2001\Toolbox\parametric_estimation.m

Listing of function C:\BrainStorm_2001\Toolbox\aic_mdl.m

function [dim, aic, mdl] = aic_mdl(s,n);
%AIC_MDL - Calculate Information Theoretic Criteria for data dimensionality
% function [dim, aic, mdl] = aic_mdl(s,n);
% Let s be the ordered eigenvalues of the data correlation matrix, calculated from
%  n samples of data.
% E.G.: F  is 151 x 375, then s = svd(F).^2; n = 375;
% Alternatively, call as [dim, aic, mdl] = aic_mdl(F);
%  where F is the spatio-temporal data matrix.
% Returns the subspace dimension, calculated by two different methods,
%  one based on Akaike (AIC), the other by Schwartz and Rissanen (MDL).
% dim is 1 x 2, with the subspace dimension as found by each method
% aic and mdl are the calculated metrics
%
% See IEEE ASSP 1985, Wax and Kailath.

%<autobegin> ---------------------- 14-Jun-2004 17:09:45 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Utility - Numeric
%
% At Check-in: $Author: Mosher $  $Revision: 9 $  $Date: 6/14/04 3:37p $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 14-Jun-2004
% 
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
%    University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
%    Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
%    CNRS, Hopital de la Salpetriere, Paris, France
% 
% See BrainStorm website at http://neuroimage.usc.edu for further information.
% 
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed  under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 14-Jun-2004 17:09:45 -----------------------

% History -----------------------------
% JCM 06-Nov-2003 creation
% -------------------------------------

if nargin == 1,
   [ignore,n] = size(s);
   s = svd(s).^2;
end

s = s/s(1); % scale to unity max for numeric handling

% test case from Wax and Kailath 1985 paper
if 0
   s = [21.2359, 2.1717, 1.4279, 1.0979, 1.0544, 0.9432, 0.7324]';
   n = 100;
end

if 0,
trim = min(find(s < .001)); % eigen values too small for practicality
if ~isempty(trim),
   s = s(1:trim);
   disp(sprintf('Trimmed to %.0f singular values',length(s)));
end
end


p = length(s);
aic = zeros(p,1);
mdl = zeros(p,1);

for k = 0:(p-1),
   
   pmk = p - k;
   
   num = prod(s((k+1):end).^(1/pmk));
   den = sum(s((k+1):end))/pmk;
   
   fac = (pmk*n)*log((num/den));
   
   aic(k+1) = -2*fac + 2*k*(2*p - k);
   
   mdl(k+1) = -fac + k*(2*p - k)*log(n)/2;
   
end

[ignore,dim] = min([aic mdl]);

dim = dim -1; % dimension of retained subspace

if 0 % debug plotting
   figure(1)
   subplot(121)
   plot([aic mdl])
   legend('aic','mdl')
   subplot(122)
   semilogy(s,'*')
end