signalstat.m

含有多种ICA算法的eeglab工具箱

% signalstat()  -  Computes and plots statistical characteristics of a signal,
%                  including the data histogram, a fitted normal distribution,
%                  a normal ditribution fitted on trimmed data, a boxplot, and
%                  the QQ-diagram. The estimates value are printed in a panel and
%                  can be read as output. Optionally, a topographic map (see TOPOPLOT)
%                  can be plotted.
%                  The boxplot and the Kolmogorov-Smirnov test require the 
%                  MATLAB Statistics Toolbox.
%
% Usage:
%   >>  signalstat( data )
%   >>  signalstat( data, plotlab, dlabel, percent );
%   >>  [M,SD,sk,k,med,zlow,zhi,tM,tSD,tndx,ksh] = ...
%                 signalstat( data, plotlab, dlabel, percent, dlabel2, map, chan_locs );
%
% Inputs:
%   data        - data vector
%
% Optional inputs:
%   plotlab     - 1: default->plot  |  0: ->no plot
%   dlabel      - A label for the data ([]: default->'Voltage [礦]')
%   percent     - percentage of data to exclude for trimmed mean & SD ([]:default->5)
%                 Excluded is 'percent'/2 high % and 'percent'/2 low %
%   dlabel2     - A title label for the statistics table
%   map         - Data vector to be displayed as topographic map. If a single integer,
%                 only the corresponding electrode location is displayed
%   chan_locs   - name of an EEG electrode position file (See >> topoplot example for format).
%                 Can also be a structure (see >> help pop_editset)
%
% Outputs:
%   M,SD        - mean and standard deviation
%   sk,k        - skewness and kurtosis
%   med         - median
%   zlow,zhi    - low and high 'percent/2'-Percentile ('percent/2'/100-Quantile)
%   tM,tSD      - trimmed mean and SD, removing data<zlow and data>zhigh
%   tndx        - index of the data retained after trimming
%   ksh         - output flag of the Kolmogorov-Smirnov test at level p=0.05 
%                 0: data could be normally distributed; 1: data are not normally distributed 
%                 -1: test could not be executed 
%
% Author: Luca Finelli, CNL / Salk Institute - SCCN, 2 August 2002
%
% See also: 
%   pop_signalstat(), qqdiagram(), eeglab() 

% Copyright (C) 2002 Luca Finelli, Salk/SCCN, La Jolla, CA

% Note: 
% QQDIAGRAM IS EQUIVALENT TO PERCENTILE/PERCENTILE PLOT
% X = EEG.data(5,:); % data
% Y = randn(1, 1000); % gaussan random distribution
% figure; qqdiagram(X, Y,  2);
% figure; plot(prctile(X,2), prctile(Y,2));

% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% $Log: signalstat.m,v $
% Revision 1.22  2003/12/03 18:49:26  arno
% msg
%
% Revision 1.21  2003/12/03 18:46:46  arno
% debug if stat toolbox abset
%
% Revision 1.20  2003/12/03 18:44:19  arno
% same
%
% Revision 1.19  2003/12/03 18:42:48  arno
% debug signal proc. toolbox absent
%
% Revision 1.18  2002/11/15 01:46:54  arno
% header for web
%
% Revision 1.17  2002/11/11 19:24:05  luca
% rearranged quantiles printout, 0.000; 1.000; ...
%
% Revision 1.16  2002/08/23 22:18:21  luca
% added topoplot (EEG channel location or component map)
%
% Revision 1.15  2002/08/22 06:39:39  luca
% uses qqdiagram(), quantile()
% excluded features if stats toolbox not present
% plots hi & low quantiles in QQ and print trimmed SD
% changed axis labels, added legend
%
% Revision 1.14  2002/08/20 18:33:56  arno
% normpdf
%
% Revision 1.13  2002/08/14 21:09:08  arno
% debugging
%
% Revision 1.12  2002/08/14 21:07:23  luca
% now checks if Statistics Toolbox is present, since required
%
% Revision 1.11  2002/08/13 22:37:44  luca
% fixed bkg color
%
% Revision 1.10  2002/08/13 18:18:14  luca
% fixed header formatting for help2html
%
% Revision 1.9  2002/08/12 20:28:56  luca
% cosmetics and title variable dlabel2 added
%
% Revision 1.8  2002/08/12 20:01:08  luca
% added the Log tag
%

function [M,SD,sk,k,med,zlow,zhi,tM,tSD,tndx,ksh] = signalstat( data, plotlab, dlabel, percent, dlabel2, map, chan_locs);

M=[]; SD=[]; sk=[]; k=[]; med=[]; zlow=[]; zhi=[]; tM=[]; tSD=[]; tndx=[]; ksh=[]; 

istats=1;	
hs = help('stats');
if isempty(hs) | ~strcmp(hs(3:20),'Statistics Toolbox')
	disp('signalstat() note: the boxplot (not shown) requires the MATLAB Statistics Toolbox');
    istats=0;
end

if (nargin<8 & nargin>5) & min(size(map))~=1
		error('signalstat(): the map input must be a vector')
end

if nargin<7 & nargin>5
	disp('signalstat(): no location file for the topographic map')
	help signalstat;
	return
end

if nargin < 6
	map = [];
end

if nargin < 5
	dlabel2 = '';
end

if nargin>3 
	if isempty(percent)
		percent=5;
	end
	if any(percent > 100) | any(percent < 0)
		error('signalstat(): percent must be between 0 and 100');
	end
end

if nargin < 4
	percent = 5;
end

if (nargin < 3 | isempty(dlabel))
	dlabel='Voltage [礦]';
end
		
if nargin < 2
  plotlab=1;
end;	
	
if ~isnumeric(plotlab)
	error('signalstat(): plotlab must be numeric');
end;

if plotlab ~= 0 & plotlab ~= 1
		error('signalstat(): plotlab must be 0 or 1');
end;
if nargin < 1
	help signalstat;
	return;
end;	

if ndims(data)>2
	error('signalstat(): data must be a vector (1-dim signal)')
end

if ~isreal(data)
	error('signalstat(): data cannot be complex')
end

fprintf('signalstat(): computing statistics...\n');

% Statistical characteristics
%----------------------------
pnts=length(data);   % number of data points
rg=max(data)-min(data);

M=mean(data);        % mean
med=median(data);    % median

vr=var(data);        % variance (N-1 normalized)
SD=std(data);        % standard deviation

if istats
	sk=skewness(data,0); % skewness (third central moment divided by
                         % the cube of the standard deviation)
    k=kurtosis(data,0);  % kurtosis (fourth central  moment divided by 
                         % fourth power of the standard deviation)
else
	sk=NaN;
	k=kurt(data);
end

% Checks on skewness and kurtosis
%--------------------------------
sklab='Distribution is symmetric';
if sk>0.01
	sklab='Distribution is right-skewed';
elseif sk < -0.01
	sklab='Distribution is left-skewed';
end

klab='';
if k>0.01
	klab='Distribution is super-Gaussian'; % i.e. kurtosis bigger then Gaussian
elseif k < -0.01
	klab='Distribution is sub-Gaussian';
end

% Estimates without the highest and lowest 'percent'/2 % of data
%---------------------------------------------------------------
pc=percent/100;

zlow = quantile(data,(pc / 2));   % low  quantile
zhi  = quantile(data,1 - pc / 2); % high quantile
tndx = find((data >= zlow & data <= zhi & ~isnan(data)));

tM=mean(data(tndx)); % mean with excluded pc/2*100% of highest and lowest values
tSD=std(data(tndx)); % trimmed SD

% Selected central tendency estimator
%------------------------------------
cte=M;

% Normal fit
%-----------
if istats
	alpha=0.05;          % 1-alpha confidence interval
	[muhat,sigmahat,muci,sigmaci] = normfit(data,alpha);
end

nbins=max(50,round(pnts/100));
[nel,binpos]=hist(data,nbins);
dx=binpos(2)-binpos(1);               % bin width

datafit=normpdf(binpos,cte,SD);       % estimated pdf