Home > chronux_1_1 > spikesort > helper > isiQuality.m

isiQuality

PURPOSE ^

ISIQUALITY Computes statistical measures of refactory period quality.

SYNOPSIS ^

function [allow, scores, cdfs] = isiQuality(unit1times, unit2times, tmin, tmax, tref, Fs)

DESCRIPTION ^

 ISIQUALITY   Computes statistical measures of refactory period quality.  
     [allow, scores, cdfs] = isiQuality(unit1times, unit2times, tmin, tmax, tref, Fs)

 Returns a boolean specifying whether the refractory period statistics are
   statistically worsened by combining the two input lists, along with two
   outputs that give more details:

 INPUTS:
   unit1times, unit2times  : (sec) original spike time lists
   tmin                    : (sec) minimum possible interval
   tmax                    : (sec) max interval for statistical comparison 
   tref                    : (sec) estimate of the refractory period
   Fs                      : (Hz)  data sampling frequency

 OUTPUTS:
   allow                   : (0 or 1) is it ok to combine the two lists?
   scores                  : (3 x 1) isi scores for [unit1 unit2 combined] lists;
                             describe the normalized fraction of spikes violating
                             the putative refractory period (1.0 means that the
                             mean density of intervals below tref is the same as the
                             mean density from tref to tmax.
   cdfs                    : (3 x m) isi cdfs for unit1, unit2 and combined lists,
                             covering the range from 0 to tmax with bin size
                             1/Fs. (i.e., m is floor(tmax*Fs)).

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [allow, scores, cdfs] = isiQuality(unit1times, unit2times, tmin, tmax, tref, Fs)
0002 
0003 % ISIQUALITY   Computes statistical measures of refactory period quality.
0004 %     [allow, scores, cdfs] = isiQuality(unit1times, unit2times, tmin, tmax, tref, Fs)
0005 %
0006 % Returns a boolean specifying whether the refractory period statistics are
0007 %   statistically worsened by combining the two input lists, along with two
0008 %   outputs that give more details:
0009 %
0010 % INPUTS:
0011 %   unit1times, unit2times  : (sec) original spike time lists
0012 %   tmin                    : (sec) minimum possible interval
0013 %   tmax                    : (sec) max interval for statistical comparison
0014 %   tref                    : (sec) estimate of the refractory period
0015 %   Fs                      : (Hz)  data sampling frequency
0016 %
0017 % OUTPUTS:
0018 %   allow                   : (0 or 1) is it ok to combine the two lists?
0019 %   scores                  : (3 x 1) isi scores for [unit1 unit2 combined] lists;
0020 %                             describe the normalized fraction of spikes violating
0021 %                             the putative refractory period (1.0 means that the
0022 %                             mean density of intervals below tref is the same as the
0023 %                             mean density from tref to tmax.
0024 %   cdfs                    : (3 x m) isi cdfs for unit1, unit2 and combined lists,
0025 %                             covering the range from 0 to tmax with bin size
0026 %                             1/Fs. (i.e., m is floor(tmax*Fs)).
0027 
0028 %   Last Modified By: sbm on Thu Jun  2 17:22:36 2005
0029 
0030 % Useful numbers
0031 f = (tref - tmin) ./ (tmax - tmin);        % fraction of range that is below tref
0032 bins = linspace(0, tmax, floor(tmax*Fs));  % using tmax*Fs bins means no info is lost to rounding
0033 refractory_bin = sum((bins <= tref));      % index of refractory period
0034 
0035 % Convert from times to interspike intervals
0036 isi1 = diff(sort(unit1times));
0037 isi2 = diff(sort(unit2times));
0038 isiT = diff(sort([unit1times; unit2times]));
0039 
0040 % Calculate histograms for the intervals below tmax.
0041 isi1hist = hist(isi1(isi1 < tmax), bins);
0042 isi2hist = hist(isi2(isi2 < tmax), bins);
0043 isiThist = hist(isiT(isiT < tmax), bins);
0044 
0045 % Count total # intervals below tmax (set 0 counts to 1 as a courtesy to the next step)
0046 isi1count = max(sum(isi1hist),1);
0047 isi2count = max(sum(isi2hist),1);
0048 isiTcount = max(sum(isiThist),1);
0049 
0050 % Make cdfs from histograms
0051 cdfs = zeros(3, length(bins));
0052 cdfs(1,:) = cumsum(isi1hist) ./ isi1count;
0053 cdfs(2,:) = cumsum(isi2hist) ./ isi2count;
0054 cdfs(3,:) = cumsum(isiThist) ./ isiTcount;
0055 
0056 %
0057 % if ((isi1count == 1) || (isi2count == 1)),  % if either original list had no violations
0058 %
0059 %
0060 % else
0061     % Compute the (scaled) difference between the initial cdfs and the combined cdfs.
0062     diffs = zeros(2, length(bins));
0063     diffs(1,:) = sqrt((isi1count*isiTcount)./(isi1count+isiTcount)) .* (cdfs(3,:) - cdfs(1,:));
0064     diffs(2,:) = sqrt((isi2count*isiTcount)./(isi2count+isiTcount)) .* (cdfs(3,:) - cdfs(2,:));
0065     
0066     % Max value of this difference in the region shorter than the refractory period
0067 
0068     dstats(1) = max(diffs(1, 1:refractory_bin));
0069     dstats(2) = max(diffs(2, 1:refractory_bin));
0070     
0071     % based on Kolmogorov-Smirnoff statistics (see Fee, Mitra, Kleinfeld, 1996)
0072     lambda = 0:0.01:2;
0073     pdfDstat = (0.5 * (1 + erf(lambda./(sqrt(2*f*(1-f)))))) - ...
0074         (0.5 .* exp(-2.*lambda.^2) .* (1 - erf((1-2*f).*lambda./(sqrt(2*f*(1-f))))));
0075     cutoff = lambda(max(find(pdfDstat < 0.95)));
0076 
0077     % Does either d-statistic exceed the statistical cutoff?
0078     if (any(dstats > cutoff)),  allow = 0;
0079     else,                       allow = 1;
0080     end
0081 % end
0082 
0083 % The cdf at the tref bin is the fraction of spikes below tref; to get a score,
0084 % we normalize by dividing by f from above.
0085 scores =  (cdfs(:, refractory_bin)) ./ f;

Generated on Sun 13-Aug-2006 11:49:44 by m2html © 2003