nt_statmatrix

stats=nt_statMatrix(x,plot_params) - calculate statistics arrays for each dim of matrix

  stats: array of statistics arrays

  x: data to stat
  plot_params: parameters

0001 function stats=nt_statMatrix(x,plot_params)
0002 %stats=nt_statMatrix(x,plot_params) - calculate statistics arrays for each dim of matrix
0003 %
0004 %  stats: array of statistics arrays
0005 %
0006 %  x: data to stat
0007 %  plot_params: parameters
0008 
0009 %% arguments
0010 if nargin<1; error('!'); end
0011 if nargin<2;
0012     plot_params.bottom=0.1;
0013     plot_params.height=0.8;
0014 end
0015 
0016 %% count true number of dimensions
0017 nDims=ndims(x);
0018 sizeX=size(x);
0019 if nDims==2 && (sizeX(1)==1 || sizeX(2)==1)
0020     nDims=1;
0021 end
0022 % if we're plotting, concatenate dimensions beyond 4th
0023 if nargout==0 && ndims(x)>4; x=x(:,:,:,:); nDims=4; end
0024 
0025 
0026 %% statistics on data
0027 stats={}; 
0028 if 0
0029     % computationally expensive
0030     for iDim=1:nDims
0031         sz=size(x);
0032         x=reshape(x,sz(1),[]);
0033         stats{iDim}.iDim=iDim;
0034         stats{iDim}.n=size(x,2);
0035         stats{iDim}.min=min(x')';
0036         stats{iDim}.max=max(x')';
0037         stats{iDim}.mean=nanmean(x')';
0038         stats{iDim}.var=nanvar(x')';
0039         x=reshape(x,sz);
0040         x=shiftdim(x,1);
0041     end
0042 else
0043     % cheaper
0044     sz=size(x);
0045     for iDim=1:nDims
0046         stats{iDim}.iDim=iDim;
0047         stats{iDim}.n=size(x,iDim);
0048         
0049         % min
0050         y=x;
0051         for k=1:iDim-1;
0052              y=squeeze(min(y,[],1));
0053         end
0054         if iDim<nDims; y=min(y(:,:),[],2); else y=min(y(:,:),[],1); end
0055         stats{iDim}.min=y(:); 
0056 
0057         % max
0058         y=x;
0059         for k=1:iDim-1;
0060             y=squeeze(max(y,[],1));
0061         end
0062         if iDim<nDims; y=max(y(:,:),[],2); else y=max(y(:,:),[],1); end
0063         stats{iDim}.max=y(:); 
0064 
0065         % mean
0066         y=x;
0067         for k=1:iDim-1;
0068             y=squeeze(nanmean(y,1));
0069         end
0070         if iDim<nDims; y=nanmean(y(:,:),2); else y=nanmean(y(:,:),1); end
0071         stats{iDim}.mean=y(:); 
0072       
0073         % var
0074         y=x.^2;
0075         for k=1:iDim-1;
0076             y=squeeze(nanmean(y,1));
0077         end
0078         if iDim<nDims; y=nanmean(y(:,:),2); else y=nanmean(y(:,:),1); end
0079         stats{iDim}.var=y(:)-stats{iDim}.mean.^2; 
0080         %nt_whoss
0081         
0082     end
0083 end
0084 
0085 
0086 %% No arguments: put up a window and plot the statistics
0087 if nargout==0;
0088     if isempty(get(0,'currentfigure'));
0089         figH=figure;
0090         set(gcf,'color',[1 1 1]);
0091     end
0092     figH=figure(gcf); %clf
0093     %set(gcf,'name',inputname(1));
0094 
0095     fontsize=12;
0096     
0097     % plot in 1,2,3 or 4 panels depending on number of dimensions
0098     switch nDims
0099         case 1
0100             axes('position',[0.05, plot_params.bottom, 0.93, plot_params.height]);  set(gca,'box','on','fontsize',fontsize);
0101             plot(x, 'k');  title(['n = ', num2str(sizeX(1))]);
0102             xlabel('samples');
0103         case 2
0104             axes('position',[0.05, plot_params.bottom, 0.45, plot_params.height]); set(gca,'box','on','fontsize',fontsize);
0105             nt_plotstats(stats{1});
0106             xlabel('samples'); title(['n = ', num2str(sizeX(1))]);
0107             axes('position',[0.53, plot_params.bottom, 0.45, plot_params.height]); set(gca,'box','on','fontsize',fontsize);
0108             nt_plotstats(stats{2});
0109             xlabel('samples'); title(['n = ', num2str(sizeX(2))]); set(gca,'ytick',[]);
0110         case 3
0111             axes('position',[0.05, plot_params.bottom, 0.3, plot_params.height]); set(gca,'box','on','fontsize',fontsize);
0112             nt_plotstats(stats{1});
0113             xlabel('samples'); title(['n = ', num2str(sizeX(1))]);
0114             axes('position',[0.37, plot_params.bottom, 0.3, plot_params.height]); set(gca,'box','on','fontsize',fontsize);
0115             nt_plotstats(stats{2});
0116             xlabel('samples'); title(['n = ', num2str(sizeX(2))]); set(gca,'ytick',[]);
0117             axes('position',[0.69, plot_params.bottom, 0.3, plot_params.height]); set(gca,'box','on','fontsize',fontsize);
0118             nt_plotstats(stats{3});
0119             xlabel('samples'); title(['n = ', num2str(sizeX(3))]); set(gca,'ytick',[]);
0120         otherwise % limit to 4 panels (last dims are concatenated)
0121             axes('position',[0.05, plot_params.bottom, 0.2, plot_params.height]); set(gca,'box','on','fontsize',fontsize);
0122             nt_plotstats(stats{1});
0123             xlabel('samples'); title(['n = ', num2str(sizeX(1))]);
0124             axes('position',[0.27, plot_params.bottom, 0.2, plot_params.height]); set(gca,'box','on','fontsize',fontsize);
0125             nt_plotstats(stats{2});
0126             xlabel('samples'); title(['n = ', num2str(sizeX(2))]); set(gca,'ytick',[]);
0127             axes('position',[0.49, plot_params.bottom, 0.2, plot_params.height]); set(gca,'box','on','fontsize',fontsize);
0128             nt_plotstats(stats{3});
0129             xlabel('samples'); title(['n = ', num2str(sizeX(3))]); set(gca,'ytick',[]);
0130             axes('position',[0.71, plot_params.bottom, 0.2, plot_params.height]); set(gca,'box','on','fontsize',fontsize);
0131             nt_plotstats(stats{4});
0132             xlabel('samples'); 
0133             title_string=['n = ', num2str(sizeX(4))];
0134             for k=5:numel(sizeX)
0135                 title_string=[title_string,'*',num2str(sizeX(k))];
0136             end
0137             title(title_string); 
0138             set(gca,'ytick',[]);
0139             
0140     end
0141 end
0142             
0143 
0144 %%
0145 function nt_plotstats(stats)
0146 %h=nt_plotstats(stats)
0147 %
0148 %  stats: stats stucture
0149 
0150 holdStatus=ishold;
0151 hold on
0152 if isfield(stats,'min')
0153     nt_plot2(1:size(stats.min,1), [stats.min,stats.max], [1 1 1]*0.9);
0154 end
0155 if isfield(stats,'var')
0156     sd=sqrt(stats.var);
0157     nt_plot2(1:size(stats.mean,1), [stats.mean+sd,stats.mean-sd], [1 1 1]* 0.5);
0158 end
0159 nt_plot2(1:size(stats.mean,1),[stats.mean,stats.mean], [1 0 0]);
0160 stats.mean(find(stats.min~=stats.max))=nan;
0161 plot(1:size(stats.mean,1),[stats.mean,stats.mean], '.b');
0162 if holdStatus;
0163     hold on;
0164 else
0165     hold off
0166 end
0167 
0168 %%
0169 function h=nt_plot2(x,y,c)
0170 %h=nt_plot2(x,y,c) - color region between two plots
0171 %
0172 %
0173 %  x: abscissa
0174 %  y: ordinate (1 or 2 columns)
0175 %  c: color (see 'fill')
0176 %
0177 %  h: vector of handles to patches
0178 %
0179 %  nt_plot2(x
0180 
0181 if nargin<1; error('!'); end
0182 
0183 % process parameters
0184 if nargin==1;
0185     y=x;
0186     x=(1:size(y,1))';
0187     c=[0 0 0];
0188 elseif nargin==2;
0189     c=[0,0,0];
0190 elseif nargin==3;
0191     ;
0192 else
0193     error('!');
0194 end
0195 
0196 % format data
0197 if size(y,2)==1
0198     y=[y,zeros(size(y))];
0199 elseif size(y,2)>2
0200     error('!');
0201 else
0202     ;
0203 end
0204 x=x(:);
0205 
0206 if 0
0207 % make sure that y(:,1)<y(:,2);
0208     yy=y;
0209     yy(:,1)=min(y,[],2);
0210     yy(:,2)=max(y,[],2);
0211     y=yy;
0212 
0213     % downsample if data array is too large
0214     TARGET_N=2000;
0215     if size(x,1)>TARGET_N
0216         DSR=ceil(size(x,1)/TARGET_N);
0217         n=floor(size(x,1)/DSR);
0218         x=x(1:DSR:end);
0219         x=x(1:n);
0220         y_extra=y(DSR*n:end,:);
0221         y=y(1:DSR*n,:);
0222         a=min(reshape(y(:,1),DSR,size(y,1)/DSR));
0223         b=max(reshape(y(:,2),DSR,size(y,1)/DSR));
0224         a(end)=min(a(end),min(y_extra(:,1)));
0225         b(end)=max(b(end),max(y_extra(:,2)));
0226         y=[a',b'];
0227     end
0228 end
0229 
0230 % draw plot
0231 yy=flipud(y(:,2));
0232 yy=yy+0.000001*max(abs(y(:))); % workaround for fill bug
0233 h=fill([x;flipud(x)],[y(:,1);yy],c, 'edgecolor', c);
0234 if x(end)>x(1); 
0235     xlim([x(1)-1,x(end)+1]);
0236 end
0237 set(gca,'box','on');
0238 
0239 %% from stats toolbox
0240 function m = nanmean(x,dim)
0241 %%
0242 
0243 % Find NaNs and set them to zero
0244 nans = isnan(x);
0245 x(nans) = 0;
0246 
0247 if nargin == 1 % let sum deal with figuring out which dimension to use
0248     % Count up non-NaNs.
0249     n = sum(~nans);
0250     n(n==0) = NaN; % prevent divideByZero warnings
0251     % Sum up non-NaNs, and divide by the number of non-NaNs.
0252     m = sum(x) ./ n;
0253 else
0254     % Count up non-NaNs.
0255     n = sum(~nans,dim);
0256     n(n==0) = NaN; % prevent divideByZero warnings
0257     % Sum up non-NaNs, and divide by the number of non-NaNs.
0258     m = sum(x,dim) ./ n;
0259 end
0260