-----------------------------------------------------------------------------
-- |
-- Module : Numeric.Statistics.Moment
-- Copyright : (c) Matthew Donadio 2002
-- License : GPL
--
-- Maintainer : m.p.donadio@ieee.org
-- Stability : experimental
-- Portability : portable
--
-- Simple module for computing the various moments of a list
--
-- Reference: Ross, NRiC
--
-----------------------------------------------------------------------------
module Numeric.Statistics.Moment (mean, var,
stddev, avgdev,
skew, kurtosis) where
-- TODO: does mean pass though the list twice? once to compute the sum,
-- and the second to compute the length?
-- TODO: does var passes through the list twice, once to compute the mean of
-- the squares, and the other to compute the mean?
-- * Functions
-- | Compute the mean of a list
--
-- @Mean(X) = 1\/N sum(i=1..N) x_i @
-- We need to use Prelude.sum intead of sum because of a buglet in the
-- Data.List library that effects nhc98
mean :: (Fractional a) => [a] -> a
mean x = Prelude.sum x / (fromIntegral.length) x
-- | Compute the variance of a list
--
-- @Var(X) = sigma^2@
--
-- @ = 1\/N-1 sum(i=1..N) (x_i-mu)^2 @
-- This is an approximation
-- var x = (mean $ map (^2) x) - mu^2
-- where mu = mean x
var :: (Fractional a) => [a] -> a
var xs = Prelude.sum (map (\x -> (x - mu)^(2::Int)) xs) / (n - 1)
where mu = mean xs
n = fromIntegral $ length $ xs
-- | Compute the standard deviation of a list
--
-- @ StdDev(X) = sigma = sqrt (Var(X)) @
stddev :: (RealFloat a) => [a] -> a
stddev x = sqrt $ var x
-- | Compute the average deviation of a list
--
-- @ AvgDev(X) = 1\/N sum(i=1..N) |x_i-mu| @
avgdev :: (RealFloat a) => [a] -> a
avgdev xs = Prelude.sum (map (\x -> abs (x - mu)) xs) / n
where mu = mean xs
n = fromIntegral $ length $ xs
-- | Compute the skew of a list
--
-- @ Skew(X) = 1\/N sum(i=1..N) ((x_i-mu)\/sigma)^3 @
skew :: (RealFloat a) => [a] -> a
skew xs = Prelude.sum (map (\x -> ((x - mu) / sigma)^(3::Int)) xs) / n
where mu = mean xs
sigma = stddev xs
n = fromIntegral $ length $ xs
-- | Compute the kurtosis of a list
--
-- @ Kurt(X) = ( 1\/N sum(i=1..N) ((x_i-mu)\/sigma)^4 ) - 3@
kurtosis :: (RealFloat a) => [a] -> a
kurtosis xs = Prelude.sum (map (\x -> ((x - mu) / sigma)^(4::Int)) xs) / n - 3
where mu = mean xs
sigma = stddev xs
n = fromIntegral $ length $ xs