{-# LANGUAGE BangPatterns, DeriveDataTypeable, DeriveGeneric,
    FlexibleContexts #-}
-- |
-- Module    : Statistics.Sample.Powers
-- Copyright : (c) 2009, 2010 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Very fast statistics over simple powers of a sample.  These can all
-- be computed efficiently in just a single pass over a sample, with
-- that pass subject to stream fusion.
--
-- The tradeoff is that some of these functions are less numerically
-- robust than their counterparts in the 'Statistics.Sample' module.
-- Where this is the case, the alternatives are noted.

module Statistics.Sample.Powers
    (
    -- * Types
      Powers

    -- * Constructor
    , powers

    -- * Descriptive functions
    , order
    , count
    , sum

    -- * Statistics of location
    , mean

    -- * Statistics of dispersion
    , variance
    , stdDev
    , varianceUnbiased

    -- * Functions over central moments
    , centralMoment
    , skewness
    , kurtosis

    -- * References
    -- $references
    ) where

import Control.Monad.ST
import Data.Aeson            (FromJSON, ToJSON)
import Data.Binary           (Binary(..))
import Data.Data             (Data, Typeable)
import Data.Vector.Binary    ()
import Data.Vector.Unboxed   ((!))
import GHC.Generics          (Generic)
import Numeric.SpecFunctions (choose)
import Prelude hiding (sum)
import Statistics.Function   (indexed)
import qualified Data.Vector          as V
import qualified Data.Vector.Generic  as G
import qualified Data.Vector.Storable as SV
import qualified Data.Vector.Unboxed  as U
import qualified Data.Vector.Unboxed.Mutable as MU
import qualified Statistics.Sample.Internal  as S

newtype Powers = Powers (U.Vector Double)
    deriving (Powers -> Powers -> Bool
(Powers -> Powers -> Bool)
-> (Powers -> Powers -> Bool) -> Eq Powers
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Powers -> Powers -> Bool
$c/= :: Powers -> Powers -> Bool
== :: Powers -> Powers -> Bool
$c== :: Powers -> Powers -> Bool
Eq, ReadPrec [Powers]
ReadPrec Powers
Int -> ReadS Powers
ReadS [Powers]
(Int -> ReadS Powers)
-> ReadS [Powers]
-> ReadPrec Powers
-> ReadPrec [Powers]
-> Read Powers
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [Powers]
$creadListPrec :: ReadPrec [Powers]
readPrec :: ReadPrec Powers
$creadPrec :: ReadPrec Powers
readList :: ReadS [Powers]
$creadList :: ReadS [Powers]
readsPrec :: Int -> ReadS Powers
$creadsPrec :: Int -> ReadS Powers
Read, Int -> Powers -> ShowS
[Powers] -> ShowS
Powers -> String
(Int -> Powers -> ShowS)
-> (Powers -> String) -> ([Powers] -> ShowS) -> Show Powers
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Powers] -> ShowS
$cshowList :: [Powers] -> ShowS
show :: Powers -> String
$cshow :: Powers -> String
showsPrec :: Int -> Powers -> ShowS
$cshowsPrec :: Int -> Powers -> ShowS
Show, Typeable, Typeable Powers
DataType
Constr
Typeable Powers
-> (forall (c :: * -> *).
    (forall d b. Data d => c (d -> b) -> d -> c b)
    -> (forall g. g -> c g) -> Powers -> c Powers)
-> (forall (c :: * -> *).
    (forall b r. Data b => c (b -> r) -> c r)
    -> (forall r. r -> c r) -> Constr -> c Powers)
-> (Powers -> Constr)
-> (Powers -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
    Typeable t =>
    (forall d. Data d => c (t d)) -> Maybe (c Powers))
-> (forall (t :: * -> * -> *) (c :: * -> *).
    Typeable t =>
    (forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c Powers))
-> ((forall b. Data b => b -> b) -> Powers -> Powers)
-> (forall r r'.
    (r -> r' -> r)
    -> r -> (forall d. Data d => d -> r') -> Powers -> r)
-> (forall r r'.
    (r' -> r -> r)
    -> r -> (forall d. Data d => d -> r') -> Powers -> r)
-> (forall u. (forall d. Data d => d -> u) -> Powers -> [u])
-> (forall u. Int -> (forall d. Data d => d -> u) -> Powers -> u)
-> (forall (m :: * -> *).
    Monad m =>
    (forall d. Data d => d -> m d) -> Powers -> m Powers)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> Powers -> m Powers)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> Powers -> m Powers)
-> Data Powers
Powers -> DataType
Powers -> Constr
(forall b. Data b => b -> b) -> Powers -> Powers
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> Powers -> c Powers
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c Powers
forall a.
Typeable a
-> (forall (c :: * -> *).
    (forall d b. Data d => c (d -> b) -> d -> c b)
    -> (forall g. g -> c g) -> a -> c a)
-> (forall (c :: * -> *).
    (forall b r. Data b => c (b -> r) -> c r)
    -> (forall r. r -> c r) -> Constr -> c a)
-> (a -> Constr)
-> (a -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
    Typeable t =>
    (forall d. Data d => c (t d)) -> Maybe (c a))
-> (forall (t :: * -> * -> *) (c :: * -> *).
    Typeable t =>
    (forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c a))
-> ((forall b. Data b => b -> b) -> a -> a)
-> (forall r r'.
    (r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall r r'.
    (r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall u. (forall d. Data d => d -> u) -> a -> [u])
-> (forall u. Int -> (forall d. Data d => d -> u) -> a -> u)
-> (forall (m :: * -> *).
    Monad m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> Data a
forall u. Int -> (forall d. Data d => d -> u) -> Powers -> u
forall u. (forall d. Data d => d -> u) -> Powers -> [u]
forall r r'.
(r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> Powers -> r
forall r r'.
(r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> Powers -> r
forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> Powers -> m Powers
forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> Powers -> m Powers
forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c Powers
forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> Powers -> c Powers
forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c Powers)
forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c Powers)
$cPowers :: Constr
$tPowers :: DataType
gmapMo :: (forall d. Data d => d -> m d) -> Powers -> m Powers
$cgmapMo :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> Powers -> m Powers
gmapMp :: (forall d. Data d => d -> m d) -> Powers -> m Powers
$cgmapMp :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> Powers -> m Powers
gmapM :: (forall d. Data d => d -> m d) -> Powers -> m Powers
$cgmapM :: forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> Powers -> m Powers
gmapQi :: Int -> (forall d. Data d => d -> u) -> Powers -> u
$cgmapQi :: forall u. Int -> (forall d. Data d => d -> u) -> Powers -> u
gmapQ :: (forall d. Data d => d -> u) -> Powers -> [u]
$cgmapQ :: forall u. (forall d. Data d => d -> u) -> Powers -> [u]
gmapQr :: (r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> Powers -> r
$cgmapQr :: forall r r'.
(r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> Powers -> r
gmapQl :: (r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> Powers -> r
$cgmapQl :: forall r r'.
(r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> Powers -> r
gmapT :: (forall b. Data b => b -> b) -> Powers -> Powers
$cgmapT :: (forall b. Data b => b -> b) -> Powers -> Powers
dataCast2 :: (forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c Powers)
$cdataCast2 :: forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c Powers)
dataCast1 :: (forall d. Data d => c (t d)) -> Maybe (c Powers)
$cdataCast1 :: forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c Powers)
dataTypeOf :: Powers -> DataType
$cdataTypeOf :: Powers -> DataType
toConstr :: Powers -> Constr
$ctoConstr :: Powers -> Constr
gunfold :: (forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c Powers
$cgunfold :: forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c Powers
gfoldl :: (forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> Powers -> c Powers
$cgfoldl :: forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> Powers -> c Powers
$cp1Data :: Typeable Powers
Data, (forall x. Powers -> Rep Powers x)
-> (forall x. Rep Powers x -> Powers) -> Generic Powers
forall x. Rep Powers x -> Powers
forall x. Powers -> Rep Powers x
forall a.
(forall x. a -> Rep a x) -> (forall x. Rep a x -> a) -> Generic a
$cto :: forall x. Rep Powers x -> Powers
$cfrom :: forall x. Powers -> Rep Powers x
Generic)

instance FromJSON Powers
instance ToJSON Powers

instance Binary Powers where
    put :: Powers -> Put
put (Powers Vector Double
v) = Vector Double -> Put
forall t. Binary t => t -> Put
put Vector Double
v
    get :: Get Powers
get = (Vector Double -> Powers) -> Get (Vector Double) -> Get Powers
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap Vector Double -> Powers
Powers Get (Vector Double)
forall t. Binary t => Get t
get

-- | O(/n/) Collect the /n/ simple powers of a sample.
--
-- Functions computed over a sample's simple powers require at least a
-- certain number (or /order/) of powers to be collected.
--
-- * To compute the /k/th 'centralMoment', at least /k/ simple powers
--   must be collected.
--
-- * For the 'variance', at least 2 simple powers are needed.
--
-- * For 'skewness', we need at least 3 simple powers.
--
-- * For 'kurtosis', at least 4 simple powers are required.
--
-- This function is subject to stream fusion.
powers :: G.Vector v Double =>
          Int                   -- ^ /n/, the number of powers, where /n/ >= 2.
       -> v Double
       -> Powers
powers :: Int -> v Double -> Powers
powers Int
k v Double
sample
  | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
2     = String -> Powers
forall a. HasCallStack => String -> a
error String
"Statistics.Sample.powers: too few powers"
  | Bool
otherwise = (forall s. ST s Powers) -> Powers
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s Powers) -> Powers)
-> (forall s. ST s Powers) -> Powers
forall a b. (a -> b) -> a -> b
$ do
      MVector s Double
acc <- Int -> Double -> ST s (MVector (PrimState (ST s)) Double)
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> a -> m (MVector (PrimState m) a)
MU.replicate Int
l Double
0
      v Double -> (Double -> ST s ()) -> ST s ()
forall (m :: * -> *) (v :: * -> *) a b.
(Monad m, Vector v a) =>
v a -> (a -> m b) -> m ()
G.forM_ v Double
sample ((Double -> ST s ()) -> ST s ()) -> (Double -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Double
x ->
        let loop :: Int -> Double -> ST s ()
loop !Int
i !Double
xk
              | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
l    = () -> ST s ()
forall (m :: * -> *) a. Monad m => a -> m a
return ()
              | Bool
otherwise = do MVector (PrimState (ST s)) Double -> Int -> Double -> ST s ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
MU.write MVector s Double
MVector (PrimState (ST s)) Double
acc Int
i (Double -> ST s ()) -> (Double -> Double) -> Double -> ST s ()
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
xk) (Double -> ST s ()) -> ST s Double -> ST s ()
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< MVector (PrimState (ST s)) Double -> Int -> ST s Double
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
MU.read MVector s Double
MVector (PrimState (ST s)) Double
acc Int
i
                               Int -> Double -> ST s ()
loop (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (Double
xk Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x)
        in Int -> Double -> ST s ()
loop Int
0 Double
1
      (Vector Double -> Powers) -> ST s (Vector Double) -> ST s Powers
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap Vector Double -> Powers
Powers (ST s (Vector Double) -> ST s Powers)
-> ST s (Vector Double) -> ST s Powers
forall a b. (a -> b) -> a -> b
$ MVector (PrimState (ST s)) Double -> ST s (Vector Double)
forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
MVector (PrimState m) a -> m (Vector a)
U.unsafeFreeze MVector s Double
MVector (PrimState (ST s)) Double
acc
  where
    l :: Int
l = Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
{-# SPECIALIZE powers :: Int -> U.Vector  Double -> Powers #-}
{-# SPECIALIZE powers :: Int -> V.Vector  Double -> Powers #-}
{-# SPECIALIZE powers :: Int -> SV.Vector Double -> Powers #-}

-- | The order (number) of simple powers collected from a 'sample'.
order :: Powers -> Int
order :: Powers -> Int
order (Powers Vector Double
pa) = Vector Double -> Int
forall a. Unbox a => Vector a -> Int
U.length Vector Double
pa Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1

-- | Compute the /k/th central moment of a sample.  The central
-- moment is also known as the moment about the mean.
centralMoment :: Int -> Powers -> Double
centralMoment :: Int -> Powers -> Double
centralMoment Int
k p :: Powers
p@(Powers Vector Double
pa)
    | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Powers -> Int
order Powers
p =
                  String -> Double
forall a. HasCallStack => String -> a
error (String
"Statistics.Sample.Powers.centralMoment: "
                         String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"invalid argument")
    | Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0    = Double
1
    | Bool
otherwise = (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
n) (Double -> Double)
-> (Vector Double -> Double) -> Vector Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Double -> Double
forall (v :: * -> *). Vector v Double => v Double -> Double
S.sum (Vector Double -> Double)
-> (Vector Double -> Vector Double) -> Vector Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Int, Double) -> Double) -> Vector (Int, Double) -> Vector Double
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
U.map (Int, Double) -> Double
go (Vector (Int, Double) -> Vector Double)
-> (Vector Double -> Vector (Int, Double))
-> Vector Double
-> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Double -> Vector (Int, Double)
forall (v :: * -> *) e.
(Vector v e, Vector v Int, Vector v (Int, e)) =>
v e -> v (Int, e)
indexed (Vector Double -> Vector (Int, Double))
-> (Vector Double -> Vector Double)
-> Vector Double
-> Vector (Int, Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Vector Double -> Vector Double
forall a. Unbox a => Int -> Vector a -> Vector a
U.take (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (Vector Double -> Double) -> Vector Double -> Double
forall a b. (a -> b) -> a -> b
$ Vector Double
pa
  where
    go :: (Int, Double) -> Double
go (Int
i , Double
e) = (Int
k Int -> Int -> Double
`choose` Int
i) Double -> Double -> Double
forall a. Num a => a -> a -> a
* ((-Double
m) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
i)) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
e
    n :: Double
n = Vector Double -> Double
forall a. Unbox a => Vector a -> a
U.head Vector Double
pa
    m :: Double
m = Powers -> Double
mean Powers
p

-- | Maximum likelihood estimate of a sample's variance.  Also known
-- as the population variance, where the denominator is /n/.  This is
-- the second central moment of the sample.
--
-- This is less numerically robust than the variance function in the
-- 'Statistics.Sample' module, but the number is essentially free to
-- compute if you have already collected a sample's simple powers.
--
-- Requires 'Powers' with 'order' at least 2.
variance :: Powers -> Double
variance :: Powers -> Double
variance = Int -> Powers -> Double
centralMoment Int
2

-- | Standard deviation.  This is simply the square root of the
-- maximum likelihood estimate of the variance.
stdDev :: Powers -> Double
stdDev :: Powers -> Double
stdDev = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> (Powers -> Double) -> Powers -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Powers -> Double
variance

-- | Unbiased estimate of a sample's variance.  Also known as the
-- sample variance, where the denominator is /n/-1.
--
-- Requires 'Powers' with 'order' at least 2.
varianceUnbiased :: Powers -> Double
varianceUnbiased :: Powers -> Double
varianceUnbiased p :: Powers
p@(Powers Vector Double
pa)
    | Double
n Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1     = Powers -> Double
variance Powers
p Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
n Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
nDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
1)
    | Bool
otherwise = Double
0
  where n :: Double
n = Vector Double -> Double
forall a. Unbox a => Vector a -> a
U.head Vector Double
pa

-- | Compute the skewness of a sample. This is a measure of the
-- asymmetry of its distribution.
--
-- A sample with negative skew is said to be /left-skewed/.  Most of
-- its mass is on the right of the distribution, with the tail on the
-- left.
--
-- > skewness . powers 3 $ U.to [1,100,101,102,103]
-- > ==> -1.497681449918257
--
-- A sample with positive skew is said to be /right-skewed/.
--
-- > skewness . powers 3 $ U.to [1,2,3,4,100]
-- > ==> 1.4975367033335198
--
-- A sample's skewness is not defined if its 'variance' is zero.
--
-- Requires 'Powers' with 'order' at least 3.
skewness :: Powers -> Double
skewness :: Powers -> Double
skewness Powers
p = Int -> Powers -> Double
centralMoment Int
3 Powers
p Double -> Double -> Double
forall a. Num a => a -> a -> a
* Powers -> Double
variance Powers
p Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (-Double
1.5)

-- | Compute the excess kurtosis of a sample.  This is a measure of
-- the \"peakedness\" of its distribution.  A high kurtosis indicates
-- that the sample's variance is due more to infrequent severe
-- deviations than to frequent modest deviations.
--
-- A sample's excess kurtosis is not defined if its 'variance' is
-- zero.
--
-- Requires 'Powers' with 'order' at least 4.
kurtosis :: Powers -> Double
kurtosis :: Powers -> Double
kurtosis Powers
p = Int -> Powers -> Double
centralMoment Int
4 Powers
p Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3
    where v :: Double
v = Powers -> Double
variance Powers
p

-- | The number of elements in the original 'Sample'.  This is the
-- sample's zeroth simple power.
count :: Powers -> Int
count :: Powers -> Int
count (Powers Vector Double
pa) = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double -> Int) -> Double -> Int
forall a b. (a -> b) -> a -> b
$ Vector Double -> Double
forall a. Unbox a => Vector a -> a
U.head Vector Double
pa

-- | The sum of elements in the original 'Sample'.  This is the
-- sample's first simple power.
sum :: Powers -> Double
sum :: Powers -> Double
sum (Powers Vector Double
pa) = Vector Double
pa Vector Double -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
! Int
1

-- | The arithmetic mean of elements in the original 'Sample'.
--
-- This is less numerically robust than the mean function in the
-- 'Statistics.Sample' module, but the number is essentially free to
-- compute if you have already collected a sample's simple powers.
mean :: Powers -> Double
mean :: Powers -> Double
mean p :: Powers
p@(Powers Vector Double
pa)
    | Double
n Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0    = Double
0
    | Bool
otherwise = Powers -> Double
sum Powers
p Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
n
    where n :: Double
n     = Vector Double -> Double
forall a. Unbox a => Vector a -> a
U.head Vector Double
pa

-- $references
--
-- * Besset, D.H. (2000) Elements of statistics.
--   /Object-oriented implementation of numerical methods/
--   ch. 9, pp. 311&#8211;331.
--   <http://www.elsevier.com/wps/product/cws_home/677916>
--
-- * Anderson, G. (2009) Compute /k/th central moments in one
--   pass. /quantblog/. <http://quantblog.wordpress.com/2009/02/07/compute-kth-central-moments-in-one-pass/>