kurtosis2(x, dimension)

function k = kurtosis2(x, dimension)
% % This program calculates the kurtosis.  
% % 
% % The normal distribution has a kurtosis of 3.
% % 
% % Description
% % 
% % If A = M x N matix, kurtosis(A) = 1 x N vector.
% % If A = M x N matix, kurtosis(A,1) = 1 x N vector.
% % If A = M x N matix, kurtosis(A,2) = M x 1 vector.
% % 
% %  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Example='1';
%
% x=randn(1,1000);      % x is the gaussian distribution 
% k = kurtosis2(x);  % returns a value close to 3
% 
% Example='2';
% 
% k = kurtosis2(x, dimension);  % returns a 1x1000 array
% dimension=1;
% 
% %  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% % Output Variables
% % 
% % k kurtosis unitless
% % 
% %  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% %     Written by William Murphy              ~2001
% %
% %     modified by Edward Zechmann 19 December 2007 
% %
% %                                 added comments    
% %                                 removed fourth moment
% %
% %     modified by Edward Zechmann 27 December 2007 
% %                                 changed filename to kurtosis2.m
% %                            
% %     modified by Edward Zechmann 29 September 2009 
% %
% %  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% % Please feel free to modify this code.
% % 

if nargin < 2
    dimension = find(size(x) > 1, 1 );
end

xsize = size(x);
if max(max(max(xsize))) > 0
    
    m = mean(x,dimension); 
    s = std(x,0,dimension);
    
    if dimension == 1
        k = (sum((x - repmat(m,[xsize(dimension),1])).^4,dimension)./(xsize(dimension)*s.^4));
    else
        %Here, we have to flip the repmat function since the dimension is different.
        k = (sum((x - repmat(m,[1,xsize(dimension)])).^4,dimension)./(xsize(dimension)*s.^4));
    end
    
else
    k = zeros(size(x));
end