# Function to calculate the unbiased population estimate or the biased sample statistic of kurtosis:

kurtosis=function(x, biased=F, na.rm=T)
{
 if (na.rm==T) x = x[!is.na(x)]
 n = length(x)
 if (n < 4)
 {
   if (na.rm==T)
   {
      cat('valid cases = ',n,'\nkurtosis is not defined for less than 4 valid cases!\n')
   }
   else
   {
      cat('cases = ',n,'\nkurtosis is not defined for less than 4 cases!\n')
   }
 }
 else
 {
   if (biased==T)
   {
      z = sqrt(n/(n-1))*scale(x)
      k = mean(z^4)-3
   }
   else
   {
      z = scale(x)
      k = sum(z^4)*n*(n+1)/((n-1)*(n-2)*(n-3))-3*(n-1)^2/((n-2)*(n-3))
   }
 k
 }
}