Accumulate - Accumulate repeated observations.
Accumulate repeated observations Y(i) = f(X1(i),X2(i),...,Xn(i)), where
the Xj(i) are binned variables (i.e. they run from 1 to Nj). This
is useful for instance to compute:
* a PSTH, where Y(i) = 1 and X(i) lists the binned spike timestamps
relative to the synchronizing events
* a JPSTH, where Y(i) = 1, X1(i) lists the binned spike timestamps of
the first neuron, and X2(i) lists those of the second neuron
* the dwell time of the animal in the environment, where Y(i) is the
sampling interval, and X1(i) and X2(i) list the binned x and y
coordinates of the animal for each video sample
* the spatial spike count of a place cell, where Y(i) = 1, and X1(i) and
X2(i) list the binned x and y coordinates of the animal for each spike
* the firing map of a place cell, as the ratio of the spatial spike count
and the dwell time
* the angular firing curve of a head direction cell, as the ratio of the
angular spike count and the dwell time
* etc.
USAGE
output = Accumulate(variables,values,<options>)
variables M x N matrix, listing the M instances of the N variables
X1(i)...Xn(i)
values optional M instances of Y(i) (circular values (e.g. angles)
must be provided in complex exponential form exp(i*theta),
where theta is in radians) (default = 1)
<options> optional list of property-value pairs (see table below)
=========================================================================
Properties Values
-------------------------------------------------------------------------
'size' output matrix size (default = max(variables))
'mode' how to deal with repeated observations: 'sum' (default),
'mean', 'min' or 'max'
=========================================================================
OUTPUT
In 'mean' mode, this function can also compute the std and 95% confidence
intervals for one or two-dimensional data (relevant only for linear data).
[mean,std,conf] = Accumulate(variables,values,'mode','mean')
In 'min' or 'max' modes, this function can also return the indices of the
min or max values:
[min,indices] = Accumulate(variables,values,'mode','min')
Missing values are set to 0 ('mean' or 'sum' modes) or NaN ('min' or 'max'
modes).
EXAMPLE
X = Bin(x,[0 1],nBins); % bin x coordinate
Y = Bin(y,[0 1],nBins); % bin y coordinate
time = Accumulate([X Y],dt); % compute occupancy map
(This is only a simplified example; to actually compute an occupancy map,
use FiringMap instead)
NOTES
1) Angular data should be provided in complex exponential form
2) Standard error and confidence intervals are valid only for linear data
SEE
See also Bin.0001 function [a,b,c] = Accumulate(variables,values,varargin) 0002 0003 %Accumulate - Accumulate repeated observations. 0004 % 0005 % Accumulate repeated observations Y(i) = f(X1(i),X2(i),...,Xn(i)), where 0006 % the Xj(i) are binned variables (i.e. they run from 1 to Nj). This 0007 % is useful for instance to compute: 0008 % 0009 % * a PSTH, where Y(i) = 1 and X(i) lists the binned spike timestamps 0010 % relative to the synchronizing events 0011 % * a JPSTH, where Y(i) = 1, X1(i) lists the binned spike timestamps of 0012 % the first neuron, and X2(i) lists those of the second neuron 0013 % * the dwell time of the animal in the environment, where Y(i) is the 0014 % sampling interval, and X1(i) and X2(i) list the binned x and y 0015 % coordinates of the animal for each video sample 0016 % * the spatial spike count of a place cell, where Y(i) = 1, and X1(i) and 0017 % X2(i) list the binned x and y coordinates of the animal for each spike 0018 % * the firing map of a place cell, as the ratio of the spatial spike count 0019 % and the dwell time 0020 % * the angular firing curve of a head direction cell, as the ratio of the 0021 % angular spike count and the dwell time 0022 % * etc. 0023 % 0024 % USAGE 0025 % 0026 % output = Accumulate(variables,values,<options>) 0027 % 0028 % variables M x N matrix, listing the M instances of the N variables 0029 % X1(i)...Xn(i) 0030 % values optional M instances of Y(i) (circular values (e.g. angles) 0031 % must be provided in complex exponential form exp(i*theta), 0032 % where theta is in radians) (default = 1) 0033 % <options> optional list of property-value pairs (see table below) 0034 % 0035 % ========================================================================= 0036 % Properties Values 0037 % ------------------------------------------------------------------------- 0038 % 'size' output matrix size (default = max(variables)) 0039 % 'mode' how to deal with repeated observations: 'sum' (default), 0040 % 'mean', 'min' or 'max' 0041 % ========================================================================= 0042 % 0043 % OUTPUT 0044 % 0045 % In 'mean' mode, this function can also compute the std and 95% confidence 0046 % intervals for one or two-dimensional data (relevant only for linear data). 0047 % 0048 % [mean,std,conf] = Accumulate(variables,values,'mode','mean') 0049 % 0050 % In 'min' or 'max' modes, this function can also return the indices of the 0051 % min or max values: 0052 % 0053 % [min,indices] = Accumulate(variables,values,'mode','min') 0054 % 0055 % Missing values are set to 0 ('mean' or 'sum' modes) or NaN ('min' or 'max' 0056 % modes). 0057 % 0058 % EXAMPLE 0059 % 0060 % X = Bin(x,[0 1],nBins); % bin x coordinate 0061 % Y = Bin(y,[0 1],nBins); % bin y coordinate 0062 % time = Accumulate([X Y],dt); % compute occupancy map 0063 % 0064 % (This is only a simplified example; to actually compute an occupancy map, 0065 % use FiringMap instead) 0066 % 0067 % NOTES 0068 % 0069 % 1) Angular data should be provided in complex exponential form 0070 % 2) Standard error and confidence intervals are valid only for linear data 0071 % 0072 % SEE 0073 % 0074 % See also Bin. 0075 0076 % Copyright (C) 2004-2016 by Michaël Zugaro, 2016 by Ralitsa Todorova, 2004 by Ken Harris 0077 % 0078 % This program is free software; you can redistribute it and/or modify 0079 % it under the terms of the GNU General Public License as published by 0080 % the Free Software Foundation; either version 3 of the License, or 0081 % (at your option) any later version. 0082 0083 % Default values 0084 outputSize = max(variables); 0085 mode = 'sum'; 0086 0087 a = []; 0088 b = []; 0089 c = []; 0090 0091 if nargin < 1, 0092 error('Incorrect number of parameters (type ''help <a href="matlab:help Accumulate">Accumulate</a>'' for details).'); 0093 end 0094 0095 if nargin < 2 0096 values = 1; 0097 end 0098 0099 % For backward compatibility: Accumulate(variables,values,outputSize) 0100 if nargin >= 3 && ~ischar(varargin{1}), 0101 outputSize = varargin{1}; 0102 varargin = {varargin{2:end}}; 0103 end 0104 0105 if isempty(variables), 0106 if ~isempty(outputSize), 0107 a = zeros(outputSize,1); 0108 end 0109 return; 0110 end 0111 0112 % Parse parameter list 0113 for i = 1:2:length(varargin), 0114 if ~ischar(varargin{i}), 0115 error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help Accumulate">Accumulate</a>'' for details).']); 0116 end 0117 switch(lower(varargin{i})), 0118 case 'size', 0119 outputSize = varargin{i+1}; 0120 if ~isivector(outputSize,'>0'), 0121 error('Incorrect value for property ''size'' (type ''help <a href="matlab:help Accumulate">Accumulate</a>'' for details).'); 0122 end 0123 case 'mode', 0124 mode = varargin{i+1}; 0125 if ~isastring(mode,'sum','mean','min','max'), 0126 error('Incorrect value for property ''mode'' (type ''help <a href="matlab:help Accumulate">Accumulate</a>'' for details).'); 0127 end 0128 otherwise, 0129 error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help Accumulate">Accumulate</a>'' for details).']); 0130 end 0131 end 0132 0133 0134 % Make sure 'values' is a Mx1 vector (possibly a scalar) 0135 if ~isdvector(values), 0136 if islvector(values), 0137 values = double(values); 0138 else 0139 error('Incorrect values (type ''help <a href="matlab:help Accumulate">Accumulate</a>'' for details)'); 0140 end 0141 end 0142 if size(values,1) == 1, values = values'; end 0143 0144 % Make sure 'variables' is a MxN matrix (possibly a Mx1 vector) of positive integers 0145 if ~isivector(variables,'>=0') & ~isimatrix(variables,'>=0'), 0146 error('Incorrect variables (type ''help <a href="matlab:help Accumulate">Accumulate</a>'' for details)'); 0147 end 0148 if size(variables,1) == 1, variables = variables'; end 0149 0150 % If 'values' is a scalar, make it a vector the same length as 'variables' 0151 if length(values) == 1, 0152 values = repmat(values,size(variables,1),1); 0153 end 0154 0155 % If 'variables' is a vector, make it a matrix 0156 if size(variables,2) == 1, 0157 variables = [variables(:) ones(length(variables),1)]; 0158 if length(outputSize) == 1, 0159 outputSize = [outputSize 1]; 0160 end 0161 end 0162 0163 % 'variables' and 'values' must have the same length 0164 if size(values,1) ~= size(variables,1), 0165 error('Incompatible sizes for variables and values (type ''help <a href="matlab:help Accumulate">Accumulate</a>'' for details)'); 0166 end 0167 0168 % Drop NaN and Inf values 0169 i = any(isnan(variables)|isinf(variables),2)|isnan(values)|isinf(values); 0170 variables(i,:) = []; 0171 values(i,:) = []; 0172 0173 % No variable should be outside the output size range 0174 if any(max(variables,[],1)>outputSize), 0175 error('Variables are out of bounds (type ''help <a href="matlab:help Accumulate">Accumulate</a>'' for details)'); 0176 end 0177 0178 % Check output size 0179 if ~isdvector(outputSize,'>0') || length(outputSize) ~= size(variables,2), 0180 error('Incorrect output size (type ''help <a href="matlab:help Accumulate">Accumulate</a>'' for details)'); 0181 end 0182 0183 % Make linear index from subscripts 0184 for i = 1:size(variables,2), 0185 subscript{i} = variables(:,i); 0186 end 0187 linearIndex = sub2ind(outputSize,subscript{:}); 0188 0189 switch mode, 0190 0191 case 'sum', 0192 tmp = sparse(linearIndex,1,values,prod(outputSize),1); 0193 a = reshape(full(tmp),outputSize); 0194 0195 case 'mean', 0196 % Divide by cell count 0197 tmp = sparse(linearIndex,1,values,prod(outputSize),1); 0198 a = reshape(full(tmp),outputSize); 0199 tmp = sparse(linearIndex,1,1,prod(outputSize),1); 0200 n = reshape(full(tmp),outputSize); 0201 a = a ./ n; 0202 a(isnan(a)) = 0; 0203 if nargout >= 2, 0204 % Standard deviation: use V = E(X2)-E(X)2 0205 tmp = sparse(linearIndex,1,values.^2,prod(outputSize),1); 0206 s = reshape(full(tmp),outputSize)./n; % E(X2) 0207 b = sqrt((s-a.^2)./(n-1).*n); % use the unbiased formula 0208 b(isnan(b)) = 0; 0209 end 0210 if nargout >= 3, 0211 % 95% confidence intervals, not computed for > 2 dimensions (too long) 0212 if length(size(variables)) <= 2, 0213 if isvector(a), 0214 for i = 1:length(a), 0215 ok = variables(:,1) == i; 0216 c(i,1) = prctile(values(ok),5); 0217 c(i,2) = prctile(values(ok),95); 0218 end 0219 else 0220 for i = 1:size(a,1), 0221 for j = 1:size(a,2), 0222 ok = variables(:,1) == i & variables(:,2) == j; 0223 c(i,j,1) = prctile(values(ok),5); 0224 c(i,j,2) = prctile(values(ok),95); 0225 end 0226 end 0227 end 0228 end 0229 end 0230 0231 case 'min', 0232 % When a vector element gets assigned multiple times, the last assignment overrides all previous assignments, 0233 % e.g. a([1 1 2]) = [4 3 6] yields a = [3 6] because a(1)=3 overrides a(1)=4. Thus, we sort values in descending 0234 % order before assigning them, so that the last one (the min) will override all others 0235 [tmp,indices] = sortrows([linearIndex values],[1 -2]); 0236 linearIndex = tmp(:,1); 0237 values = tmp(:,2); 0238 a = nan([outputSize,1]); 0239 a(linearIndex) = values; 0240 a = reshape(a,outputSize); 0241 if nargout >= 2, 0242 b = nan(outputSize); 0243 b(linearIndex) = indices; 0244 end 0245 0246 case 'max', 0247 % When a vector element gets assigned multiple times, the last assignment overrides all previous assignments, 0248 % e.g. a([1 1 2]) = [4 3 6] yields a = [3 6] because a(1)=3 overrides a(1)=4. Thus, we sort values in ascending 0249 % order before assigning them, so that the last one (the max) will override all others 0250 [tmp,indices] = sortrows([linearIndex values],[1 2]); 0251 linearIndex = tmp(:,1); 0252 values = tmp(:,2); 0253 a = nan([outputSize,1]); 0254 a(linearIndex) = values; 0255 a = reshape(a,outputSize); 0256 if nargout >= 2, 0257 b = nan(outputSize); 0258 b(linearIndex) = indices; 0259 end 0260 0261 end 0262 0263