function [rcum] = cumulants(data,cumulant_start,cumulant_end)
    n = 1;
    cum(n) = 0;
    n = 2;
    while(n <= cumulant_end)        
        kn = moment(data,n);
        k = 1;        
        while(k <= n-1)                       
            nCk = prod(k:(n-1))/prod(1:(n-k)); %actually it's 'n-1 C k-1'
            kn = kn-nCk*cum(1,k)*moment(data,(n-k));
            k = k+1;           
        end
        cum(n) = kn;
        n = n+1;            
    end       
    rcum = cum(cumulant_start:1:cumulant_end);
end