MapStats - Compute statistics for a map Z = f(X,Y), or a curve Z = f(X).
Compute statistics on a continuous map, where one time-varying variable Z
is represented as a function of one or two time-varying variables X and Y.
The variable Z can either be a point process (typically, a list of spike
timestamps) or a continuous measure (e.g. the instantaneous velocity of
the animal, the spectral power of an LFP channel in a given frequency band,
the coherence between two oscillating LFP channels, etc.) Typical examples
of X and Y include spatial coordinates and angular directions.
USAGE
stats = MapStats(map,<options>)
map map obtained using <a href="matlab:help Map">Map</a>
<options> optional list of property-value pairs (see table below)
=========================================================================
Properties Values
-------------------------------------------------------------------------
'threshold' values above threshold*peak belong to the field
(default = 0.2)
'minSize' fields smaller than this size are considered spurious
and ignored (default = 100 for 2D, 10 for 1D)
'minPeak' peaks smaller than this size are considered spurious
and ignored (default = 1)
'type' 'll' if X and Y are linear, 'cl' if X is circular and Y
linear, 'lc' if X is linear and Y circular, or 'cc' if X
and Y are circular - for 1D data, a single letter is used
(default = 'll')
'verbose' display processing information (default = 'off')
=========================================================================
OUTPUT
stats.x abscissa of the maximum value (in bins)
stats.y ordinate of the maximum value (in bins)
stats.peak in-field maximum value
stats.mean in-field mean value
stats.size field size (in bins)
stats.field field (1 = bin in field, 0 = bin not in field)
stats.fieldX field x boundaries (in bins)
stats.fieldY field y boundaries (in bins)
stats.specificity spatial specificity (in bits, see Skaggs et al., 1993)
For 1D circular data:
stats.m mean angle
stats.mode distribution mode (in bins)
stats.r mean resultant length
stats.k von Mises concentration
SEE
See also Map, FiringMap, FiringCurve.0001 function stats = MapStats(map,varargin) 0002 0003 %MapStats - Compute statistics for a map Z = f(X,Y), or a curve Z = f(X). 0004 % 0005 % Compute statistics on a continuous map, where one time-varying variable Z 0006 % is represented as a function of one or two time-varying variables X and Y. 0007 % The variable Z can either be a point process (typically, a list of spike 0008 % timestamps) or a continuous measure (e.g. the instantaneous velocity of 0009 % the animal, the spectral power of an LFP channel in a given frequency band, 0010 % the coherence between two oscillating LFP channels, etc.) Typical examples 0011 % of X and Y include spatial coordinates and angular directions. 0012 % 0013 % USAGE 0014 % 0015 % stats = MapStats(map,<options>) 0016 % 0017 % map map obtained using <a href="matlab:help Map">Map</a> 0018 % <options> optional list of property-value pairs (see table below) 0019 % 0020 % ========================================================================= 0021 % Properties Values 0022 % ------------------------------------------------------------------------- 0023 % 'threshold' values above threshold*peak belong to the field 0024 % (default = 0.2) 0025 % 'minSize' fields smaller than this size are considered spurious 0026 % and ignored (default = 100 for 2D, 10 for 1D) 0027 % 'minPeak' peaks smaller than this size are considered spurious 0028 % and ignored (default = 1) 0029 % 'type' 'll' if X and Y are linear, 'cl' if X is circular and Y 0030 % linear, 'lc' if X is linear and Y circular, or 'cc' if X 0031 % and Y are circular - for 1D data, a single letter is used 0032 % (default = 'll') 0033 % 'verbose' display processing information (default = 'off') 0034 % ========================================================================= 0035 % 0036 % OUTPUT 0037 % 0038 % stats.x abscissa of the maximum value (in bins) 0039 % stats.y ordinate of the maximum value (in bins) 0040 % stats.peak in-field maximum value 0041 % stats.mean in-field mean value 0042 % stats.size field size (in bins) 0043 % stats.field field (1 = bin in field, 0 = bin not in field) 0044 % stats.fieldX field x boundaries (in bins) 0045 % stats.fieldY field y boundaries (in bins) 0046 % stats.specificity spatial specificity (in bits, see Skaggs et al., 1993) 0047 % 0048 % For 1D circular data: 0049 % 0050 % stats.m mean angle 0051 % stats.mode distribution mode (in bins) 0052 % stats.r mean resultant length 0053 % stats.k von Mises concentration 0054 % 0055 % SEE 0056 % 0057 % See also Map, FiringMap, FiringCurve. 0058 0059 % Copyright (C) 2002-2012 by Michaël Zugaro 0060 % 0061 % This program is free software; you can redistribute it and/or modify 0062 % it under the terms of the GNU General Public License as published by 0063 % the Free Software Foundation; either version 3 of the License, or 0064 % (at your option) any later version. 0065 0066 % Check number of parameters 0067 if nargin < 1 | mod(length(varargin),2) ~= 0, 0068 error('Incorrect number of parameters (type ''help <a href="matlab:help MapStats">MapStats</a>'' for details).'); 0069 end 0070 0071 if ~isfield(map,'z'), 0072 map.z = map.rate; 0073 map = rmfield(map,'rate'); 0074 end 0075 0076 % Default values 0077 threshold = 0.2; 0078 minPeak = 1; 0079 type = 'll'; 0080 verbose = 0; 0081 0082 nDims = sum(size(map.z)>=2); 0083 if nDims == 2, 0084 minSize = 100; 0085 else 0086 minSize = 10; 0087 end 0088 0089 % Parse parameter list 0090 for i = 1:2:length(varargin), 0091 if ~ischar(varargin{i}), 0092 error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help MapStats">MapStats</a>'' for details).']); 0093 end 0094 switch(lower(varargin{i})), 0095 0096 case 'threshold', 0097 threshold = varargin{i+1}; 0098 if ~isdscalar(threshold,'>=0','<=1'), 0099 error('Incorrect value for property ''threshold'' (type ''help <a href="matlab:help MapStats">MapStats</a>'' for details).'); 0100 end 0101 0102 case 'minsize', 0103 minSize = varargin{i+1}; 0104 if ~isdscalar(minSize,'>=0'), 0105 error('Incorrect value for property ''minSize'' (type ''help <a href="matlab:help MapStats">MapStats</a>'' for details).'); 0106 end 0107 0108 case 'minpeak', 0109 minPeak = varargin{i+1}; 0110 if ~isdscalar(minPeak,'>=0'), 0111 error('Incorrect value for property ''minPeak'' (type ''help <a href="matlab:help MapStats">MapStats</a>'' for details).'); 0112 end 0113 0114 case 'type', 0115 type = lower(varargin{i+1}); 0116 if (nDims == 1 && ~isastring(type,'c','l')) || (nDims == 2 && ~isastring(type,'cc','cl','lc','ll')), 0117 error('Incorrect value for property ''type'' (type ''help <a href="matlab:help MapStats">MapStats</a>'' for details).'); 0118 end 0119 0120 case {'verbose','debug'}, 0121 verbose = lower(varargin{i+1}); 0122 if ~isastring(verbose,'on','off'), 0123 error('Incorrect value for property ''verbose'' (type ''help <a href="matlab:help MapStats">MapStats</a>'' for details).'); 0124 end 0125 verbose = strcmp(verbose,'on'); 0126 0127 otherwise, 0128 error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help MapStats">MapStats</a>'' for details).']); 0129 0130 end 0131 end 0132 0133 % Are X and/or Y circular? 0134 circX = size(map.z,2) > 1 && strcmp(type(1),'c'); 0135 circY = size(map.z,1) > 1 && ((size(map.z,2) > 1 && strcmp(type(2),'c')) || strcmp(type(1),'c')); 0136 0137 % Default values 0138 stats.x = NaN; 0139 stats.y = NaN; 0140 stats.field = logical(zeros(0,0,0)); 0141 stats.size = 0; 0142 stats.peak = 0; 0143 stats.mean = 0; 0144 stats.fieldX = [NaN NaN]; 0145 stats.fieldY = [NaN NaN]; 0146 stats.specificity = 0; 0147 stats.m = nan; 0148 stats.r = nan; 0149 stats.mode = nan; 0150 stats.k = nan; 0151 0152 if isempty(map.z), return; end 0153 0154 % Compute the spatial specificity of the map, based on the formula proposed by Skaggs et al. (1993): 0155 % specificity = SUM { p(i) . lambda(i)/lambda . log2(lambda(i)/lambda) } 0156 T = sum(map.time(:)); 0157 p_i = map.time/(T+eps); 0158 lambda_i = map.z; 0159 lambda = lambda_i(:)'*p_i(:); 0160 if T == 0 || lambda == 0, 0161 stats.specificity = 0; 0162 else 0163 stats.specificity = sum(sum(p_i.*lambda_i/lambda.*log2(lambda_i/lambda))); 0164 end 0165 0166 % Determine the field as the connex area around the peak where the value or rate is > threshold*peak 0167 % There can be two or more fields 0168 0169 if max(max(map.z)) == 0, 0170 stats.field = logical(zeros(size(map.z))); 0171 return; 0172 end 0173 0174 nBinsX = max([1 length(map.x)]); % minimum number of bins is 1 0175 nBinsY = max([1 length(map.y)]); 0176 0177 % Each time we find a field, we will remove it from the map; make a copy first 0178 z = map.z; 0179 % Try to find more fields until no remaining bin exceeds min value 0180 i = 1; 0181 while true, 0182 % Are there any candidate (unvisited) peaks left? 0183 [peak,idx] = max(z(:)); 0184 if peak < minPeak, break; end 0185 % Determine coordinates of largest candidate peak 0186 [y,x] = ind2sub(size(z),idx); 0187 % Find field (using min threshold for inclusion) 0188 field1 = FindField(z,x,y,peak*threshold,circX,circY); 0189 size1 = sum(field1(:)); 0190 % Does this field include two coalescent subfields? 0191 % To answer this question, we simply re-run the same field-searching procedure on the field 0192 % we then either keep the original field or choose the subfield if the latter is less than 0193 % 1/2 the size of the former 0194 m = peak*threshold; 0195 field2 = FindField(z-m,x,y,(peak-m)*threshold,circX,circY); 0196 size2 = sum(field2(:)); 0197 if size2< 1/2*size1, 0198 field = field2; 0199 tc = ' ';sc = '*'; % for debugging messages 0200 else 0201 field = field1; 0202 tc = '*';sc = ' '; % for debugging messages 0203 end 0204 % Display debugging info 0205 if verbose, 0206 disp([int2zstr(i,2) ') peak ' num2str(peak) ' @ (' int2str(x) ',' int2str(y) ')']); 0207 disp([' ' tc ' field size ' int2str(size1)]); 0208 disp([' ' sc ' subfield size ' int2str(size2)]); 0209 disp(' '); 0210 figure; 0211 if nDims == 1, 0212 plot(z);hold on; 0213 PlotIntervals(ToIntervals(field1),'rectangles'); 0214 PlotIntervals(ToIntervals(field2),'bars'); 0215 ylabel(tc); 0216 else 0217 subplot(3,1,1);imagesc(z);xlabel('Data'); 0218 subplot(3,1,2);imagesc(field1);clim([0 1]);xlabel('Field'); 0219 subplot(3,1,3);imagesc(field2);clim([0 1]);ylabel(tc);xlabel('Subfield'); 0220 end 0221 end 0222 fieldSize = sum(field(:)); 0223 % Keep this field if its size is sufficient 0224 if fieldSize > minSize, 0225 stats.field(:,:,i) = field; 0226 stats.size(i) = fieldSize; 0227 stats.peak(i) = peak; 0228 stats.mean(i) = mean(z(field)); 0229 idx = find(field & z == peak); 0230 [stats.y(i),stats.x(i)] = ind2sub(size(z),idx(1)); 0231 [x,y] = FieldBoundaries(field,circX,circY); 0232 [stats.fieldX(i,:),stats.fieldY(i,:)] = FieldBoundaries(field,circX,circY); 0233 i = i + 1; 0234 end 0235 % Mark field bins as visited 0236 z(field) = NaN; 0237 if all(isnan(z)), break; end 0238 end 0239 0240 % ----------------------------- Circular statistics ----------------------------- 0241 0242 if ~(nDims == 1 && circX), return; end 0243 0244 complex = map.z .* exp(j*map.x*2*pi) / sum(map.z); 0245 stats.m = angle(nanmean(complex)); 0246 stats.r = abs(nansum(complex)); 0247 % [~,stats.mode] = max(map.z); 0248 [maxZ, iMax] = max(map.z); 0249 TMax = map.x(iMax); 0250 stats.mode = TMax*2*pi; 0251 0252 % Compute kappa (Fisher, "Statistical Analysis of Circular Data" p.88) 0253 if stats.r < 0.53, 0254 stats.k = 2*stats.r+stats.r^3+5*stats.r^5/6; 0255 elseif stats.r < 0.85, 0256 stats.k = -0.4+1.39*stats.r+0.43/(1-stats.r); 0257 else 0258 stats.k = 1/(stats.r^3-4*stats.r^2+3*stats.r); 0259 end 0260 0261 % ------------------------------- Helper functions ------------------------------- 0262 0263 % Field boundaries (circumscribed rectangle) 0264 0265 function [x,y] = FieldBoundaries(field,circX,circY) 0266 0267 % Find boundaries 0268 x = find(any(field,1)); 0269 if isempty(x), 0270 x = [NaN NaN]; 0271 else 0272 x = [x(1) x(end)]; 0273 end 0274 y = find(any(field,2)); 0275 if isempty(y), 0276 y = [NaN NaN]; 0277 else 0278 y = [y(1) y(end)]; 0279 end 0280 0281 % The above works in almost all cases; it fails however for circular coordinates if the field extends 0282 % around an edge, e.g. for angles between 350° and 30° 0283 0284 if circX && x(1) == 1 && x(2) == size(field,2), 0285 xx = find(~all(field,1)); 0286 if ~isempty(xx), 0287 x = [xx(end) xx(1)]; 0288 end 0289 end 0290 if circY && y(1) == 1 && y(2) == size(field,1), 0291 yy = find(~all(field,2)); 0292 if ~isempty(yy), 0293 y = [yy(end) yy(1)]; 0294 end 0295 end 0296