rollingstats

classdef rollingstats
    %RUNNINGSTATS Rolling descriptive statistics
    %   This class is a convenience utility made to compute descriptive
    %   statistics (mean, variance, kurtosis, etc...) of a data sample in a
    %   rolling fashion.
    
    properties (SetAccess = private)
        n = 0;         % Number of elements added
        mean = 0;
        M2 = 0;
        M3 = 0;
        M4 = 0;
        min = [];
        max = [];
    end
    
    %% Dependent properties
    properties (Dependent = true, SetAccess = private, Transient = true)
        variance
        std
        stderr
        kurtosis
        skewness
    end
    
    methods
        
        function v = get.variance(obj)
            v = obj.M2 / obj.n;
        end
        
        function v = get.std(obj)
            v = sqrt( double(obj.variance) );
        end
        
        function v = get.stderr(obj)
            v = obj.std / sqrt(obj.n);
        end
        
        function v = get.kurtosis(obj)
            v = (obj.n.*obj.M4) ./ (obj.M2.*obj.M2) - 3;
        end
        
        function v = get.skewness(obj)
            v = sqrt(obj.n) * obj.M3 ./ realpow(obj.M2, 1.5);
        end
    end
    
    %% Public methods
    
    methods
        
        function obj = add(obj, x)
            
            n1 = obj.n;
            obj.n = n1 + 1;
            delta = x - obj.mean;
            delta_n = delta / obj.n;
            delta_n2 = delta_n .* delta_n;
            term1 = delta .* delta_n * n1;
            obj.mean = obj.mean + delta_n;
            
            obj.M4 = obj.M4 + ...
                term1 .* delta_n2 * (obj.n*obj.n - 3*obj.n + 3) + ...
                6 * delta_n2 .* obj.M2 - ...
                4 * delta_n .* obj.M3;
            
            obj.M3 = obj.M3 + ...
                term1 .* delta_n .* (obj.n - 2) - ...
                3 * delta_n .* obj.M2;
            
            obj.M2 = obj.M2 + term1;
            
            if n1 == 0
                
                obj.min = x;
                obj.max = x;
            else
                
                
                indexmin = x < obj.min;
                obj.min(indexmin) = x(indexmin);
                
                indexmax = x > obj.max;
                obj.max(indexmax) = x(indexmax);
            end
            
            
        end
        
        
    end
    
end