SyncHist

SyncHist - Compute a histogram on event-synchronized samples (e.g. a PSTH).

 Compute a histogram on event-synchronized point-process (e.g. spikes) or
 continuous (e.g. voltages) samples, obtained using the function <a href="matlab:help Sync">Sync</a>.
 Individual values Xi(t) (for trial i and time t) can be averaged or
 summed across trials. Alternatively, complete time-varying distributions
 can also be computed, i.e. this function can compute for each time t the
 distribution of Xi(t) across trials i.

  USAGE

    To sum across trials:
    [sum,timeBins] = SyncHist(synchronized,indices,'mode','sum',<options>)

    To average across trials (computes a frequency for point processes):
    [mean,error,timeBins] = SyncHist(synchronized,indices,'mode','mean',<options>)
    (see 'error' property)

    To compute the full distribution:
    [distribution,timeBins,valueBins] = SyncHist(synchronized,indices,'mode','dist',<options>)

    synchronized   resynchronized samples (obtained using <a href="matlab:help Sync">Sync</a>)
    indices        synchronizing event indices (obtained using <a href="matlab:help Sync">Sync</a>)
    <options>      optional list of property-value pairs (see table below)

    =========================================================================
     Properties    Values
    -------------------------------------------------------------------------
     'mode'        type of histogram: 'sum' (default), 'mean' or 'dist'
     'durations'   durations before and after synchronizing events for each
                   trial (in s) (default = [-0.5 0.5])
     'nBins'       total number of time bins (default 100)
     'type'        either 'linear' or 'circular' (use radians) depending on
                   the type of data in the raster (default 'linear')
     'smooth'      standard deviation for Gaussian kernel (default 0, no
                   smoothing)
     'bins'        [m M n], lower and upper bounds, and number of bins,
                   respectively (default [min max 100]) (only for 'dist' mode)
     'error'       either 'std' to compute standard deviation (default),
                   'sem' to compute standard error of the mean, or '95%' to
                   compute 95% confidence intervals (only in 'mean' mode)
    =========================================================================

  EXAMPLE 1

    % Compute and plot spike raster
    [raster,indices] = Sync(spikes,stimuli);
    figure;PlotSync(raster,indices);

    % Compute and plot PSTH
    [m,t] = SyncHist(raster,indices,'mode','sum');
    figure;bar(t,m);

  EXAMPLE 2

    % Make velocity timestamps relative to trial onsets
    [vsync,indices] = Sync(velocity,trials);

    % Compute velocity distribution across trials
    [d,x,t] = SyncHist(vsync,indices,'mode','dist');

    % Plot distribution as color map
    figure;PlotColorMap(d,1);

  SEE

    See also Sync, SyncMap, PlotSync, PETHTransition.

0001 function [a,b,c] = SyncHist(synchronized,indices,varargin)
0002 
0003 %SyncHist - Compute a histogram on event-synchronized samples (e.g. a PSTH).
0004 %
0005 % Compute a histogram on event-synchronized point-process (e.g. spikes) or
0006 % continuous (e.g. voltages) samples, obtained using the function <a href="matlab:help Sync">Sync</a>.
0007 % Individual values Xi(t) (for trial i and time t) can be averaged or
0008 % summed across trials. Alternatively, complete time-varying distributions
0009 % can also be computed, i.e. this function can compute for each time t the
0010 % distribution of Xi(t) across trials i.
0011 %
0012 %  USAGE
0013 %
0014 %    To sum across trials:
0015 %    [sum,timeBins] = SyncHist(synchronized,indices,'mode','sum',<options>)
0016 %
0017 %    To average across trials (computes a frequency for point processes):
0018 %    [mean,error,timeBins] = SyncHist(synchronized,indices,'mode','mean',<options>)
0019 %    (see 'error' property)
0020 %
0021 %    To compute the full distribution:
0022 %    [distribution,timeBins,valueBins] = SyncHist(synchronized,indices,'mode','dist',<options>)
0023 %
0024 %    synchronized   resynchronized samples (obtained using <a href="matlab:help Sync">Sync</a>)
0025 %    indices        synchronizing event indices (obtained using <a href="matlab:help Sync">Sync</a>)
0026 %    <options>      optional list of property-value pairs (see table below)
0027 %
0028 %    =========================================================================
0029 %     Properties    Values
0030 %    -------------------------------------------------------------------------
0031 %     'mode'        type of histogram: 'sum' (default), 'mean' or 'dist'
0032 %     'durations'   durations before and after synchronizing events for each
0033 %                   trial (in s) (default = [-0.5 0.5])
0034 %     'nBins'       total number of time bins (default 100)
0035 %     'type'        either 'linear' or 'circular' (use radians) depending on
0036 %                   the type of data in the raster (default 'linear')
0037 %     'smooth'      standard deviation for Gaussian kernel (default 0, no
0038 %                   smoothing)
0039 %     'bins'        [m M n], lower and upper bounds, and number of bins,
0040 %                   respectively (default [min max 100]) (only for 'dist' mode)
0041 %     'error'       either 'std' to compute standard deviation (default),
0042 %                   'sem' to compute standard error of the mean, or '95%' to
0043 %                   compute 95% confidence intervals (only in 'mean' mode)
0044 %    =========================================================================
0045 %
0046 %  EXAMPLE 1
0047 %
0048 %    % Compute and plot spike raster
0049 %    [raster,indices] = Sync(spikes,stimuli);
0050 %    figure;PlotSync(raster,indices);
0051 %
0052 %    % Compute and plot PSTH
0053 %    [m,t] = SyncHist(raster,indices,'mode','sum');
0054 %    figure;bar(t,m);
0055 %
0056 %  EXAMPLE 2
0057 %
0058 %    % Make velocity timestamps relative to trial onsets
0059 %    [vsync,indices] = Sync(velocity,trials);
0060 %
0061 %    % Compute velocity distribution across trials
0062 %    [d,x,t] = SyncHist(vsync,indices,'mode','dist');
0063 %
0064 %    % Plot distribution as color map
0065 %    figure;PlotColorMap(d,1);
0066 %
0067 %  SEE
0068 %
0069 %    See also Sync, SyncMap, PlotSync, PETHTransition.
0070 
0071 % Copyright (C) 2004-2013 by Michaël Zugaro
0072 %
0073 % This program is free software; you can redistribute it and/or modify
0074 % it under the terms of the GNU General Public License as published by
0075 % the Free Software Foundation; either version 3 of the License, or
0076 % (at your option) any later version.
0077 
0078 % Default values
0079 durations = [-0.5 0.5];
0080 nBins = 100;
0081 smooth = 0;
0082 type = 'linear';
0083 bins = [];
0084 a = [];
0085 b = [];
0086 c = [];
0087 mode = 'sum';
0088 error = 'std';
0089 
0090 % Check number of parameters
0091 if nargin < 2 | mod(length(varargin),2) ~= 0,
0092   error('Incorrect number of parameters (type ''help <a href="matlab:help SyncHist">SyncHist</a>'' for details).');
0093 end
0094 
0095 if isempty(indices) || isempty(synchronized), return; end
0096 
0097 % Check parameter sizes
0098 if size(indices,2) ~= 1,
0099     error('Parameter ''indices'' is not a vector (type ''help <a href="matlab:help SyncHist">SyncHist</a>'' for details).');
0100 end
0101 if size(indices,1) ~= size(synchronized,1),
0102     error('Parameters ''synchronized'' and ''indices'' have different lengths (type ''help <a href="matlab:help SyncHist">SyncHist</a>'' for details).');
0103 end
0104 
0105 % Parse parameter list
0106 for i = 1:2:length(varargin),
0107     if ~ischar(varargin{i}),
0108         error(['Parameter ' num2str(i+1) ' is not a property (type ''help <a href="matlab:help SyncHist">SyncHist</a>'' for details).']);
0109     end
0110     switch(lower(varargin{i})),
0111         case 'mode',
0112             mode = lower(varargin{i+1});
0113             if ~isastring(mode,'sum','mean','dist'),
0114                 error('Incorrect value for property ''mode'' (type ''help <a href="matlab:help SyncHist">SyncHist</a>'' for details).');
0115             end
0116 
0117         case 'durations',
0118             durations = varargin{i+1};
0119             if ~isdvector(durations,'#2','<'),
0120                 error('Incorrect value for property ''durations'' (type ''help <a href="matlab:help SyncHist">SyncHist</a>'' for details).');
0121             end
0122 
0123         case 'nbins',
0124             nBins = varargin{i+1};
0125             if ~isdscalar(nBins,'>0'),
0126                 error('Incorrect value for property ''nBins'' (type ''help <a href="matlab:help SyncHist">SyncHist</a>'' for details).');
0127             end
0128 
0129         case 'smooth',
0130             smooth = varargin{i+1};
0131             if ~isdvector(smooth,'>=0') | length(smooth) > 2,
0132                 error('Incorrect value for property ''smooth'' (type ''help <a href="matlab:help SyncHist">SyncHist</a>'' for details).');
0133             end
0134 
0135         case {'bins','hist'},
0136             bins = varargin{i+1};
0137             if ~isdvector(bins,'#3') | bins(1) > bins(2) | ~isiscalar(bins(3),'>0'),
0138                 error('Incorrect value for property ''bins'' (type ''help <a href="matlab:help SyncHist">SyncHist</a>'' for details).');
0139             end
0140 
0141         case 'type',
0142             type = lower(varargin{i+1});
0143             if ~isastring(type,'linear','circular'),
0144                 error('Incorrect value for property ''type'' (type ''help <a href="matlab:help SyncHist">SyncHist</a>'' for details).');
0145             end
0146 
0147         case 'error',
0148             error = lower(varargin{i+1});
0149             if ~isastring(error,'std','sem','95%'),
0150                 error('Incorrect value for property ''error'' (type ''help <a href="matlab:help SyncHist">SyncHist</a>'' for details).');
0151             end
0152 
0153         otherwise,
0154             error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help SyncHist">SyncHist</a>'' for details).']);
0155 
0156     end
0157 end
0158 
0159 % Check number of output parameters
0160 switch(mode),
0161     case 'sum',
0162         if nargout > 2,
0163             error('Incorrect number of output parameters (type ''help <a href="matlab:help SyncHist">SyncHist</a>'' for details).');
0164         end
0165     case 'mean',
0166         if nargout > 4,
0167             error('Incorrect number of output parameters (type ''help <a href="matlab:help SyncHist">SyncHist</a>'' for details).');
0168         end
0169     case 'dist',
0170         if nargout > 3,
0171             error('Incorrect number of output parameters (type ''help <a href="matlab:help SyncHist">SyncHist</a>'' for details).');
0172         end
0173 end
0174 
0175 if size(synchronized,2) == 1,
0176     % Point process data
0177     pointProcess = true;
0178     synchronized(:,2) = 1;
0179     if isempty(bins),
0180         bins = [0 max(synchronized(:,2)) 100];
0181     end
0182 else
0183     % Continuous data
0184     pointProcess = false;
0185     if isempty(bins),
0186         bins = [min(min(synchronized(:,2:end))) max(max(synchronized(:,2:end))) 100];
0187     end
0188 end
0189 
0190 % Compute histogram
0191 nTrials = max(indices);
0192 start = durations(1);
0193 stop = durations(2);
0194 timeBinSize = (stop - start)/nBins;
0195 timeBins = (start:timeBinSize:stop-timeBinSize)+timeBinSize/2;
0196 binnedTime = Bin(synchronized(:,1),[start stop],nBins,'trim');
0197 
0198 if strcmp(type,'circular'),
0199     range = isradians(synchronized(:,2:end));
0200     for j = 2:size(synchronized,2),
0201         synchronized(:,j) = exp(i*synchronized(:,j));
0202     end
0203 end
0204 
0205 switch(mode),
0206 
0207     case 'sum',
0208         % Smoothed sum
0209         for j = 2:size(synchronized,2),
0210             a(:,j-1) = Smooth(Accumulate(binnedTime,synchronized(:,j),nBins),smooth);
0211         end
0212         b = timeBins;
0213         c = [];
0214     case 'mean',
0215         if pointProcess,
0216             % Smoothed mean
0217             for j = 2:size(synchronized,2),
0218                 a(:,j-1) = Smooth(Accumulate(binnedTime,synchronized(:,j),nBins),smooth)/(nTrials*timeBinSize);
0219             end
0220             b = []; % Not yet implemented
0221         else
0222             % Number of values in each time bin
0223             n = Smooth(Accumulate(binnedTime,1,nBins),smooth);
0224             % Smoothed mean
0225             for j = 2:size(synchronized,2),
0226                 a(:,j-1) = Smooth(Accumulate(binnedTime,synchronized(:,j),nBins),smooth)./n;
0227             end
0228             if strcmp(type,'circular'),
0229                 b = []; % Error not computed for circular data
0230             else
0231                 if strcmp(error,'std'),
0232                     % Variance
0233                     for j = 2:size(synchronized,2),
0234                         sq = Smooth(Accumulate(binnedTime,synchronized(:,j).^2,nBins),smooth)./n; % E(X²)
0235                         b(:,j-1) = sqrt(sq-a(:,j-1).^2).*n./(n-1);
0236                     end
0237                 elseif strcmp(error,'sem'),
0238                     % Standard error of the mean (also use V = E(X²)-E(X)²)
0239                     for j = 2:size(synchronized,2),
0240                         sq = Smooth(Accumulate(binnedTime,synchronized(:,j).^2,nBins),smooth)./n; % E(X²)
0241                         b(:,j-1) = sqrt(sq-a(:,j-1).^2).*n./(n-1)./sqrt(n);
0242                     end
0243                 elseif strcmp(error,'95%'),
0244                     % 95% confidence intervals
0245                     for j = 2:size(synchronized,2),
0246                         for i = 1:nBins,
0247                             b(i,j-1,1) = prctile(synchronized(binnedTime==i,2),5);
0248                             b(i,j-1,2) = prctile(synchronized(binnedTime==i,2),95);
0249                         end
0250                     end
0251                 end
0252             end
0253         end
0254         c = timeBins;
0255     case 'dist',
0256         if pointProcess,
0257             % Number of points in each (time,trial) bin
0258             n = Accumulate([binnedTime indices],1);
0259             % Smoothed distribution
0260             a = Smooth(Accumulate([n(:)+1 repmat((1:nBins)',nTrials,1)],1),smooth)/nTrials;
0261             b = timeBins;
0262             c = 0:max(max(n));
0263         else
0264             minValue = bins(1);
0265             maxValue = bins(2);
0266             nValueBins = bins(3);
0267             % Number of values in each time bin
0268             n = Accumulate(binnedTime,1,nBins);
0269             for j = 2:size(synchronized,2),
0270                 % Binned values
0271                 binnedValues = Bin(synchronized(:,j),[minValue maxValue],nValueBins);
0272                 % Smoothed distribution
0273                 a(:,:,j-1) = Smooth(Accumulate([binnedValues binnedTime],1,[nValueBins nBins]),smooth)./Smooth(repmat(n',nValueBins,1),smooth);
0274             end
0275             range = maxValue - minValue;
0276             binSize = range/nValueBins;
0277             b = timeBins;
0278             c = (minValue:range/nValueBins:maxValue-binSize)'+binSize/2;
0279         end
0280 end
0281 
0282 if strcmp(type,'circular'),
0283     a(:,2:end) = wrap(a(:,2:end),range);
0284 end