Home > NoiseTools > nt_statmatrix.m

# nt_statmatrix

## PURPOSE

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

## SYNOPSIS

function stats=nt_statMatrix(x,plot_params)

## DESCRIPTION

    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

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- nt_dataview [p,data]=nt_dataview(data,p) - view data sets

## SUBFUNCTIONS

- function nt_plotstats(stats)
- function h=nt_plot2(x,y,c)
- function m = nanmean(x,dim)

## SOURCE CODE

    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