signalstat.m

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

datafit=datafit*pnts*dx;

tdatafit=normpdf(binpos,cte,tSD);     % estimated pdf with trimmed SD
tdatafit=tdatafit*pnts*dx;

%datarnd=normrnd(cte,sigmahat,1,pnts); % synthetic data

if istats
	% Goodness-of-fit hypothesis test
	%--------------------------------
	kstail = 0; % 0 = 2-sided test
	
	CDF=normcdf(data,cte,sigmahat);        % estimated cdf
	
	[ksh,ksp,ksstat,kscv] = kstest(data,[data', CDF'],alpha,kstail); % Kolmogorov-Smirnov test
	
	kstestlab='Kolmogorov-Smirnov test: verified Gaussian';
	kscol=[0.2 1 0.2];
	if ksh
		kstestlab='Kolmogorov-Smirnov test: not Gaussian';
		kscol=[.7 .3 .3];
	end 
end

% Graphics
%-------------------------------------

if plotlab
  figure
  try, icadefs; set(gcf, 'color', BACKCOLOR); catch, end;
  COLOR = [0.56 .66 .9];
  set(gcf,'NumberTitle','off','Name','Signal statistics -- signalstat()')
  fwidth=800;  % figure size in pixels
  fheight=600;
    
  funits=get(0,'Units');
  set(0,'Units','Pixel')
  scnsize=get(0,'ScreenSize');
  fpos=[round((scnsize(3)-fwidth)/2),round((scnsize(4)-fheight)/2),fwidth,fheight];
  set(0,'Units',funits)
   
  set(gcf,'Position',fpos)
 
  % Plotting the histogram
  %------------------------
  subplot(2,2,1)
  
  hist(data,nbins);  % 'XLim',[-125 125]
  uapos = get(gca,'Position');
  xlim  = get(gca,'XLim');
  set(gca,'FontSize',14)
  xlabel(dlabel)
  title('Data Histogram and Fitted Normal PDF')
 
  % HM=pnts*normpdf(M,muhat,sigmahat)/2; % FWHM height
  % plot([cte-sd cte+sd],[HM HM],'r--','LineWidth',2)

  % Overplotting a normal distribution
  %-----------------------------------
  hold on
  h1=plot(binpos,datafit,'c','LineWidth',2);
  set(gca,'XLim',xlim)
  
  % Overplotting a normal distribution from trimmed SD
  %---------------------------------------------------
  h2=plot(binpos,tdatafit,'y');
  ymin=get(gca,'YLim');
  plot([zlow zlow],[0 ymin(2)/20],'y','LineWidth',2) % low  percentile
  plot([zhi  zhi], [0 ymin(2)/20],'y','LineWidth',2) % high percentile
  set(gca,'XLim',xlim)

  % Overplotting a mean and zero line
  %----------------------------------
 % xmin=get(gca,'XLim');
 % ymin=get(gca,'YLim');
  plot([0 0],ymin,'k')
  h3=plot([cte cte],ymin,'r--','LineWidth',2);
  set(gca,'Color',COLOR,'XMinorTick','on','XLim',xlim)
  
  if istats
	  set(gca,'XTick',[])	  
  elseif ~istats & strcmp(dlabel,'Voltage [礦]')
	  set(gca,'XTick',[-125, -75, -25, 0, 25, 75,  125],...
			  'XTickLabel',['-125' ; ' -75' ; ' -25' ; '  0 ' ; ' 25 ' ; ' 75 ' ; ' 125'])
  end
  
  if strcmp(dlabel,'Voltage [礦]')
	  set(gca,'XLim',[-125 125])
  end

  set(gca,'FontSize',10)
  H=[h1 h2 h3];
  legend(H,'Gaussian fit','Trimmed G.fit','Mean') 
  legend boxoff
  
  zoom off

  % Boxplot
  %--------------------
  if istats
	  subplot(2,2,3)
	  
	  boxplot(data,1,'+',0,1.5)
	  lapos=get(gca,'Position');
	  set(gca,'Position', [uapos(1) uapos(2)-uapos(4)/3 uapos(3) uapos(4)/3])
	  nlapos=get(gca,'Position');
	  hold on
	  ymin2=get(gca,'YLim');
	  plot([0 0],[0 ymin2(2)],'k')
	  plot([cte cte],[0 ymin(2)],'r--','LineWidth',2)
	  set(gca,'FontSize',14,'XMinorTick','on') 
      set(gca,'XLim',xlim)
	  
      if strcmp(dlabel,'Voltage [礦]')
		  set(gca,'XTick',[-125 -75 -25 0 25 75  125],...
				  'XTickLabel',['-125' ; ' -75' ; ' -25' ; '  0 ' ; ' 25 ' ; ' 75 ' ; ' 125'],...
				  'XLim',[-125 125])
	  end
	 
	  xlabel(dlabel)
	  ylabel('')
	  zoom off
  end
  
  % QQ plot
  %--------
  subplot(2,2,2)
  
  qqdiagram(data)
  apos=get(gca,'Position');
  set(gca,'Position', [1-uapos(1)-uapos(3) uapos(2)-uapos(4)/3 (uapos(3)) (uapos(4)+uapos(4)/3)])
  set(gca,'XTick',[-4 -2 0 2 4])
  xmin=get(gca,'XLim');
  ymin=get(gca,'YLim');
  hold on
  plot([xmin(1) xmin(1)+diff(xmin)/20],[zlow zlow],'y-','LineWidth',2)
  plot([xmin(1) xmin(1)+diff(xmin)/20],[zhi  zhi] ,'y-','LineWidth',2)
  set(gca,'XLim',xmin);
  %plot([0 0],ymin,'k--')
  set(gca,'FontSize',14)
  xlabel('Standard Normal Quantiles [Std.Dev.]')
  if strcmp(dlabel,'Voltage [礦]')
	  ylabel('Ordered Observations [礦]')
  elseif strcmp(dlabel,'Component Activity')
	  ylabel('Ordered Observations [rel. 礦]')
  else
	  ylabel('Ordered Observations')
  end
  
  title('QQ Plot (Data vs Standard Normal)')
  set(gca,'Color',COLOR)
  
  % TOPO plot   
  %---------
  if (~isempty(map))
	  sbplot(7,9,6)
	  % th=axes('Position',[]);
	  % subplot('Position',[.10 .86 .20 .14]); 
	  fprintf('signalstat(): plotting a topographic map...\n');
	  if length(map) == 1
		  topoplot(map,chan_locs,'electrodes','off', ...
				   'style', 'blank', 'emarkersize1chan', 10);
	  else
		  topoplot(map,chan_locs,'electrodes','off');
	  end;
	  axis('square')
  end 

  % Color schemes
  %--------------
  nero    = [0 0 0];
  bordeau = [0.7 0.3 0.3];
  rosso   = [1 0 0];
  roschi1 = [1 .3 0.4];
  giallo1 = [1 .9 0];
  arancio = [1 .5 .3];
  verchi1 = [0.2 1 0.2];
  verchi2 = [0.7 1 0.7];
  verchi3 = [0.4 1 0.4];
  verscu1 = [0.1 0.7 0.2];
  bluchi1 = [0.4 0.9 1];
  bluchi2 = [.2 .5 .7];
  grichia = [.95 .95 .95];
  
  % Data axis
  %------------
  bgcolor=COLOR;
  
  dah = axes('Position',[uapos(1) .1 1-2*uapos(1) .25]);
  set(dah,'Box','on','Color',bgcolor,'XTick',[],'YTick',[],'FontName','Courier','FontSize',12,'FontWeight','demi')
  title(dlabel2,'FontWeight','bold','FontSize',14,'FontName','Arial');
  text(0.05,0.9,['Mean:          ' num2str(M,3)]  ,'Color',bordeau,'FontName','Courier','FontSize',12,'FontWeight','demi')
  text(0.05,0.8,['Trimmed mean:  ' num2str(tM,3)] ,'Color',giallo1,'FontName','Courier','FontSize',12,'FontWeight','demi')

  text(0.05,0.5,['Standard dev.: ' num2str(SD,4)] ,'Color',bordeau,'FontName','Courier','FontSize',12,'FontWeight','demi')
  text(0.05,0.4,['Trimmed st.d.: ' num2str(tSD,4)],'Color',giallo1,'FontName','Courier','FontSize',12,'FontWeight','demi')
  text(0.05,0.3,['Variance:      ' num2str(vr,4)] ,'Color',giallo1,'FontName','Courier','FontSize',12,'FontWeight','demi')
  text(0.05,0.2,['Range:         ' num2str(rg,4)] ,'Color',giallo1,'FontName','Courier','FontSize',12,'FontWeight','demi') 
  text(0.05,0.1,['Data points:   ' num2str(pnts)] ,'Color',bordeau,'FontName','Courier','FontSize',12,'FontWeight','demi')
   
  text(0.4,0.9,[num2str(percent/2/100,'%1.3f') '-quantile: ' num2str(zlow,3)] ,'Color',verchi2,...
	   'FontName','Courier','FontSize',12,'FontWeight','demi')
  text(0.4,0.8,['0.5  -quantile: ',num2str(med,3),'  (median)'],'Color',verchi1,'FontName','Courier','FontSize',12,'FontWeight','demi')
  text(0.4,0.7,[num2str((100-percent/2)/100,'%1.3f') '-quantile:  ' num2str(zhi,3)] ,'Color',verchi2,...
	   'FontName','Courier','FontSize',12,'FontWeight','demi')
  
  text(0.4,0.3,['Kurtosis: ' num2str(k, 3) ' (0 if Gaussian)'] ,'Color',verchi1,...
				'FontName','Courier','FontSize',12,'FontWeight','demi')
  text(0.4,0.2,klab,'Color',verchi2,'FontName','Courier','FontSize',12,'FontWeight','demi')
 
  if istats
	  text(0.4,0.5,['Skewness: ' num2str(sk,3) ' (0 if Gaussian)'] ,'Color',verchi1,...
		   			'FontName','Courier','FontSize',12,'FontWeight','demi')
	  text(0.4,0.4,sklab,'Color',verchi2,'FontName','Courier','FontSize',12,'FontWeight','demi')
      text(0.4,0.1,kstestlab,'Color',kscol,'FontName','Courier','FontSize',12,'FontWeight','demi')
  end
  axcopy;
end

%--------------------------------------------------
% clone of the normpdf function of the stat toolbox
function fitvals = normpdf(myvals,mymean,mystd)
if nargin < 3,
    mystd = 1;
end
if nargin < 2;
    mymean = 0;
end
if length(mymean) < length(myvals)
	tmpmean = mymean;
	mymean = zeros(size(myvals));
	mymean(:) = tmpmean;
end;
if length(mystd) < length(myvals)
	tmpmean = mystd;
	mystd = zeros(size(myvals));
	mystd(:) = tmpmean;
end;
mymean(1:10);
mystd(1:10);

fitvals = zeros(size(myvals));
tmp = find(mystd > 0);
if any(tmp)
    myvalsn = (myvals(tmp) - mymean(tmp)) ./ mystd(tmp);
    fitvals(tmp) = exp(-0.5 * myvalsn .^2) ./ (sqrt(2*pi) .* mystd(tmp));
end
tmp1 = find(mystd <= 0);
if any(tmp1)
    tmp2   = NaN;
    fitvals(tmp1) = tmp2(ones(size(tmp1)));
end