Home > atmlab > math > running_stats.m

running_stats

PURPOSE ^

running_stats Calculate running statistics

SYNOPSIS ^

function stat_vectors = running_stats(x, y, binsize, stat_handles)

DESCRIPTION ^

 running_stats Calculate running statistics

 FORMAT

   stats = running_stat(x, y, binsize, stat_handles)

 IN

   x [vector]
       independent quantity, must be constantly spaced

   y [vector]
       quantity to calculate statistics for

   winsize [scalar]
       size of window.
       if positive, interpreted as same unit as x
       if negative but >-1, interpreted as fraction (negative is just a
       trick here to pass a flag), then extra stat_handle appended showing
       number of samples used for each statistic
       

   stat_handles [cell array of function handles]
       functions to apply to each running window
       Note that each function MUST accept a dimension as its second
       argument, such as median(X, 2)!

 OUT

   stat_vectors [cell array of vectors]
       Output of all handles in stat_handles for each running window.
       Padded with nans at begin and end where sliding stats are
       undefined.

 See also: smooth, medfilt1, movingstd, filter

 Note: smooth is in the curve-fitting toolbox, medfilt1 is in the signal
 processing toolbox.
 * in the curve fitting toolbox, otherwise see
   http://stackoverflow.com/a/1516168/974555
 ^ from File Exchange

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

running_stats.m

SOURCE CODE ^

0001 function stat_vectors = running_stats(x, y, binsize, stat_handles)
0002 % running_stats Calculate running statistics
0003 %
0004 % FORMAT
0005 %
0006 %   stats = running_stat(x, y, binsize, stat_handles)
0007 %
0008 % IN
0009 %
0010 %   x [vector]
0011 %       independent quantity, must be constantly spaced
0012 %
0013 %   y [vector]
0014 %       quantity to calculate statistics for
0015 %
0016 %   winsize [scalar]
0017 %       size of window.
0018 %       if positive, interpreted as same unit as x
0019 %       if negative but >-1, interpreted as fraction (negative is just a
0020 %       trick here to pass a flag), then extra stat_handle appended showing
0021 %       number of samples used for each statistic
0022 %
0023 %
0024 %   stat_handles [cell array of function handles]
0025 %       functions to apply to each running window
0026 %       Note that each function MUST accept a dimension as its second
0027 %       argument, such as median(X, 2)!
0028 %
0029 % OUT
0030 %
0031 %   stat_vectors [cell array of vectors]
0032 %       Output of all handles in stat_handles for each running window.
0033 %       Padded with nans at begin and end where sliding stats are
0034 %       undefined.
0035 %
0036 % See also: smooth, medfilt1, movingstd, filter
0037 %
0038 % Note: smooth is in the curve-fitting toolbox, medfilt1 is in the signal
0039 % processing toolbox.
0040 % * in the curve fitting toolbox, otherwise see
0041 %   http://stackoverflow.com/a/1516168/974555
0042 % ^ from File Exchange
0043 
0044 % $Id: running_stats.m 7901 2012-09-27 15:38:16Z gerrit $
0045 
0046 % convert binsize [x-unit] to binsize [n-index]
0047 
0048 % some input checks
0049 
0050 rqre_same_size(x, y);
0051 
0052 % verify that x-unit is constantly spaced
0053 % error otherwise
0054 
0055 df = diff(x);
0056 spacing_xunit = unique(df);
0057 if length(spacing_xunit)>1
0058     error(['atmlab:' mfilename ':notunique'], ...
0059         ['Error in calculation of running statistics: ' ...
0060         'spacing in independent quantity not constant']);
0061 end
0062 
0063 
0064 if binsize > 0
0065     stat_vectors = running_stats_linear(x, y, binsize, stat_handles, spacing_xunit);
0066 elseif binsize < 0 && binsize > -1
0067     % fractional binsize, e.g. 0.01 is all within 1%
0068     stat_vectors = running_stats_fractional(x, y, binsize, stat_handles, spacing_xunit);
0069 else
0070     error(['atmlab:' mfilename], ...
0071         'Binsize must be >0 or <0 & >-1. Found: %.3f', binsize);
0072 end
0073 
0074 end
0075 
0076 function stat_vectors = running_stats_linear(x, y, binsize_xunit, stat_handles, spacing_xunit)
0077     
0078 
0079 N = round(binsize_xunit/spacing_xunit);
0080 L = length(x);
0081 
0082 % pre-allocate
0083 stat_vectors = preallocate_stat_vectors(x, stat_handles);
0084 
0085 % Yes: idea from
0086 % http://www.mathworks.com/matlabcentral/fileexchange/29029-summarizes-data-using-a-sliding-window-without-loops/content/movingstat.m
0087 
0088 % re-arrange data so there are columns for each window (needs memory)
0089 data_arranged = y(repmat(0:N-1,L-N+1,1)+cumsum(ones(L-N+1,N)));
0090 if iseven(N)
0091     ran = (N/2):(L-N/2);
0092 else
0093     ran = (ceil(N/2)):(L-floor(N/2));
0094 end
0095 
0096 % TODO: do in chunks similar to medfilt1
0097 for j = 1:length(stat_handles)
0098     stat_vectors{j}(ran) = stat_handles{j}(data_arranged, 2);
0099 end
0100 
0101 % this is the slow version:
0102 % for i = (binsize_n/2+1):(length(x)-binsize_n/2)
0103 %     slice = (i-(binsize_n/2)):(i+(binsize_n/2));
0104 %     segment = y(slice);
0105 %     for j = 1:length(stat_handles)
0106 %         stat_vectors{j}(i) = stat_handles{j}(segment);
0107 %     end
0108 % end
0109 
0110 
0111 end
0112 
0113 function stat_vectors = running_stats_fractional(x, y, binsize_frac, stat_handles, spacing_xunit)
0114 
0115 L = length(x);
0116 
0117 % fraction is passed negative
0118 binsize_frac = -binsize_frac;
0119 
0120 % always tell the length
0121 stat_handles = [stat_handles {@size}];
0122 
0123 % pre-allocate
0124 stat_vectors = preallocate_stat_vectors(x, stat_handles);
0125 
0126 % get lower and upper edge of segments for each x
0127 
0128 lo = x .* (1 - binsize_frac/2);
0129 hi = x .* (1 + binsize_frac/2);
0130 
0131 % translate this to indices
0132 
0133 N = ceil(2 * (x * binsize_frac) / spacing_xunit);
0134 fst = find(N > 100, 1, 'first'); % arbitrary start
0135 
0136 % try brute-force first (POITROAE)
0137 for i = fst:length(y)
0138     % search for indices corresponding to lo(i) and hi(i)
0139     slice = floor(i-(N(i)/2)):ceil(i+(N(i)/2));
0140     if (slice(1) < 1) || (slice(end) > L)
0141         continue
0142     end
0143     segment = y(slice);
0144     
0145     for j = 1:length(stat_handles)
0146         stat_vectors{j}(i) = stat_handles{j}(segment, 1);
0147     end
0148     
0149     if mod(i, 100000)==0
0150         logtext(atmlab('OUT'), '%d/%d done\n', i, length(y));
0151     end
0152 end
0153 
0154 end
0155 
0156 function stat_vectors = preallocate_stat_vectors(x, stat_handles)
0157 
0158 stat_vectors = cell(size(stat_handles));
0159 for i = 1:length(stat_handles)
0160     stat_vectors{i} = nan(size(x));
0161 end
0162 
0163 end

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005