Home > chronux_2_00 > 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 % Useful numbers
0029 f = (tref - tmin) ./ (tmax - tmin);        % fraction of range that is below tref
0030 bins = linspace(0, tmax, floor(tmax*Fs));  % using tmax*Fs bins means no info is lost to rounding
0031 refractory_bin = sum((bins <= tref));      % index of refractory period
0032 
0033 % Convert from times to interspike intervals
0034 isi1 = diff(sort(unit1times));
0035 isi2 = diff(sort(unit2times));
0036 isiT = diff(sort([unit1times; unit2times]));
0037 
0038 % Calculate histograms for the intervals below tmax.
0039 isi1hist = hist(isi1(isi1 < tmax), bins);
0040 isi2hist = hist(isi2(isi2 < tmax), bins);
0041 isiThist = hist(isiT(isiT < tmax), bins);
0042 
0043 % Count total # intervals below tmax (set 0 counts to 1 as a courtesy to the next step)
0044 isi1count = max(sum(isi1hist),1);
0045 isi2count = max(sum(isi2hist),1);
0046 isiTcount = max(sum(isiThist),1);
0047 
0048 % Make cdfs from histograms
0049 cdfs = zeros(3, length(bins));
0050 cdfs(1,:) = cumsum(isi1hist) ./ isi1count;
0051 cdfs(2,:) = cumsum(isi2hist) ./ isi2count;
0052 cdfs(3,:) = cumsum(isiThist) ./ isiTcount;
0053 
0054 %
0055 % if ((isi1count == 1) || (isi2count == 1)),  % if either original list had no violations
0056 %
0057 %
0058 % else
0059     % Compute the (scaled) difference between the initial cdfs and the combined cdfs.
0060     diffs = zeros(2, length(bins));
0061     diffs(1,:) = sqrt((isi1count*isiTcount)./(isi1count+isiTcount)) .* (cdfs(3,:) - cdfs(1,:));
0062     diffs(2,:) = sqrt((isi2count*isiTcount)./(isi2count+isiTcount)) .* (cdfs(3,:) - cdfs(2,:));
0063     
0064     % Max value of this difference in the region shorter than the refractory period
0065 
0066     dstats(1) = max(diffs(1, 1:refractory_bin));
0067     dstats(2) = max(diffs(2, 1:refractory_bin));
0068     
0069     % based on Kolmogorov-Smirnoff statistics (see Fee, Mitra, Kleinfeld, 1996)
0070     lambda = 0:0.01:2;
0071     pdfDstat = (0.5 * (1 + erf(lambda./(sqrt(2*f*(1-f)))))) - ...
0072         (0.5 .* exp(-2.*lambda.^2) .* (1 - erf((1-2*f).*lambda./(sqrt(2*f*(1-f))))));
0073     cutoff = lambda(max(find(pdfDstat < 0.95)));
0074 
0075     % Does either d-statistic exceed the statistical cutoff?
0076     if (any(dstats > cutoff)),  allow = 0;
0077     else                        allow = 1;
0078     end
0079 % end
0080 
0081 % The cdf at the tref bin is the fraction of spikes below tref; to get a score,
0082 % we normalize by dividing by f from above.
0083 scores =  (cdfs(:, refractory_bin)) ./ f;

Generated on Fri 15-Aug-2008 11:35:42 by m2html © 2003