{-# LANGUAGE FlexibleContexts #-}
-- |
-- Module    : Statistics.Sample
-- Copyright : (c) 2008 Don Stewart, 2009 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Commonly used sample statistics, also known as descriptive
-- statistics.

module Statistics.Sample
    (
    -- * Types
      Sample
    , WeightedSample
    -- * Descriptive functions
    , range

    -- * Statistics of location
    , mean
    , welfordMean
    , meanWeighted
    , harmonicMean
    , geometricMean

    -- * Statistics of dispersion
    -- $variance

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

    -- ** Two-pass functions (numerically robust)
    -- $robust
    , variance
    , varianceUnbiased
    , meanVariance
    , meanVarianceUnb
    , stdDev
    , varianceWeighted
    , stdErrMean

    -- ** Single-pass functions (faster, less safe)
    -- $cancellation
    , fastVariance
    , fastVarianceUnbiased
    , fastStdDev

    -- * Joint distributions
    , covariance
    , correlation
    , pair
    -- * References
    -- $references
    ) where

import Statistics.Function (minMax)
import Statistics.Sample.Internal (robustSumVar, sum)
import Statistics.Types.Internal  (Sample,WeightedSample)
import qualified Data.Vector as V
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U

-- Operator ^ will be overridden
import Prelude hiding ((^), sum)

-- | /O(n)/ Range. The difference between the largest and smallest
-- elements of a sample.
range :: (G.Vector v Double) => v Double -> Double
range :: forall (v :: * -> *). Vector v Double => v Double -> Double
range v Double
s = Double
hi forall a. Num a => a -> a -> a
- Double
lo
    where (Double
lo , Double
hi) = forall (v :: * -> *).
Vector v Double =>
v Double -> (Double, Double)
minMax v Double
s
{-# INLINE range #-}

-- | /O(n)/ Arithmetic mean.  This uses Kahan-Babuška-Neumaier
-- summation, so is more accurate than 'welfordMean' unless the input
-- values are very large.
mean :: (G.Vector v Double) => v Double -> Double
mean :: forall (v :: * -> *). Vector v Double => v Double -> Double
mean v Double
xs = forall (v :: * -> *). Vector v Double => v Double -> Double
sum v Double
xs forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v Double
xs)
{-# SPECIALIZE mean :: U.Vector Double -> Double #-}
{-# SPECIALIZE mean :: V.Vector Double -> Double #-}

-- | /O(n)/ Arithmetic mean.  This uses Welford's algorithm to provide
-- numerical stability, using a single pass over the sample data.
--
-- Compared to 'mean', this loses a surprising amount of precision
-- unless the inputs are very large.
welfordMean :: (G.Vector v Double) => v Double -> Double
welfordMean :: forall (v :: * -> *). Vector v Double => v Double -> Double
welfordMean = T -> Double
fini forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (v :: * -> *) b a.
Vector v b =>
(a -> b -> a) -> a -> v b -> a
G.foldl' T -> Double -> T
go (Double -> Int -> T
T Double
0 Int
0)
  where
    fini :: T -> Double
fini (T Double
a Int
_) = Double
a
    go :: T -> Double -> T
go (T Double
m Int
n) Double
x = Double -> Int -> T
T Double
m' Int
n'
        where m' :: Double
m' = Double
m forall a. Num a => a -> a -> a
+ (Double
x forall a. Num a => a -> a -> a
- Double
m) forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n'
              n' :: Int
n' = Int
n forall a. Num a => a -> a -> a
+ Int
1
{-# SPECIALIZE welfordMean :: U.Vector Double -> Double #-}
{-# SPECIALIZE welfordMean :: V.Vector Double -> Double #-}

-- | /O(n)/ Arithmetic mean for weighted sample. It uses a single-pass
-- algorithm analogous to the one used by 'welfordMean'.
meanWeighted :: (G.Vector v (Double,Double)) => v (Double,Double) -> Double
meanWeighted :: forall (v :: * -> *).
Vector v (Double, Double) =>
v (Double, Double) -> Double
meanWeighted = V -> Double
fini forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (v :: * -> *) b a.
Vector v b =>
(a -> b -> a) -> a -> v b -> a
G.foldl' V -> (Double, Double) -> V
go (Double -> Double -> V
V Double
0 Double
0)
    where
      fini :: V -> Double
fini (V Double
a Double
_) = Double
a
      go :: V -> (Double, Double) -> V
go (V Double
m Double
w) (Double
x,Double
xw) = Double -> Double -> V
V Double
m' Double
w'
          where m' :: Double
m' | Double
w' forall a. Eq a => a -> a -> Bool
== Double
0   = Double
0
                   | Bool
otherwise = Double
m forall a. Num a => a -> a -> a
+ Double
xw forall a. Num a => a -> a -> a
* (Double
x forall a. Num a => a -> a -> a
- Double
m) forall a. Fractional a => a -> a -> a
/ Double
w'
                w' :: Double
w' = Double
w forall a. Num a => a -> a -> a
+ Double
xw
{-# INLINE meanWeighted #-}

-- | /O(n)/ Harmonic mean.  This algorithm performs a single pass over
-- the sample.
harmonicMean :: (G.Vector v Double) => v Double -> Double
harmonicMean :: forall (v :: * -> *). Vector v Double => v Double -> Double
harmonicMean = T -> Double
fini forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (v :: * -> *) b a.
Vector v b =>
(a -> b -> a) -> a -> v b -> a
G.foldl' T -> Double -> T
go (Double -> Int -> T
T Double
0 Int
0)
  where
    fini :: T -> Double
fini (T Double
b Int
a) = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
a forall a. Fractional a => a -> a -> a
/ Double
b
    go :: T -> Double -> T
go (T Double
x Int
y) Double
n = Double -> Int -> T
T (Double
x forall a. Num a => a -> a -> a
+ (Double
1forall a. Fractional a => a -> a -> a
/Double
n)) (Int
yforall a. Num a => a -> a -> a
+Int
1)
{-# INLINE harmonicMean #-}

-- | /O(n)/ Geometric mean of a sample containing no negative values.
geometricMean :: (G.Vector v Double) => v Double -> Double
geometricMean :: forall (v :: * -> *). Vector v Double => v Double -> Double
geometricMean = forall a. Floating a => a -> a
exp forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (v :: * -> *). Vector v Double => v Double -> Double
mean forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (v :: * -> *) a b.
(Vector v a, Vector v b) =>
(a -> b) -> v a -> v b
G.map forall a. Floating a => a -> a
log
{-# INLINE geometricMean #-}

-- | Compute the /k/th central moment of a sample.  The central moment
-- is also known as the moment about the mean.
--
-- This function performs two passes over the sample, so is not subject
-- to stream fusion.
--
-- For samples containing many values very close to the mean, this
-- function is subject to inaccuracy due to catastrophic cancellation.
centralMoment :: (G.Vector v Double) => Int -> v Double -> Double
centralMoment :: forall (v :: * -> *). Vector v Double => Int -> v Double -> Double
centralMoment Int
a v Double
xs
    | Int
a forall a. Ord a => a -> a -> Bool
< Int
0  = forall a. HasCallStack => [Char] -> a
error [Char]
"Statistics.Sample.centralMoment: negative input"
    | Int
a forall a. Eq a => a -> a -> Bool
== Int
0 = Double
1
    | Int
a forall a. Eq a => a -> a -> Bool
== Int
1 = Double
0
    | Bool
otherwise = forall (v :: * -> *). Vector v Double => v Double -> Double
sum (forall (v :: * -> *) a b.
(Vector v a, Vector v b) =>
(a -> b) -> v a -> v b
G.map Double -> Double
go v Double
xs) forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v Double
xs)
  where
    go :: Double -> Double
go Double
x = (Double
xforall a. Num a => a -> a -> a
-Double
m) Double -> Int -> Double
^ Int
a
    m :: Double
m    = forall (v :: * -> *). Vector v Double => v Double -> Double
mean v Double
xs
{-# SPECIALIZE centralMoment :: Int -> U.Vector Double -> Double #-}
{-# SPECIALIZE centralMoment :: Int -> V.Vector Double -> Double #-}

-- | Compute the /k/th and /j/th central moments of a sample.
--
-- This function performs two passes over the sample, so is not subject
-- to stream fusion.
--
-- For samples containing many values very close to the mean, this
-- function is subject to inaccuracy due to catastrophic cancellation.
centralMoments :: (G.Vector v Double) => Int -> Int -> v Double -> (Double, Double)
centralMoments :: forall (v :: * -> *).
Vector v Double =>
Int -> Int -> v Double -> (Double, Double)
centralMoments Int
a Int
b v Double
xs
    | Int
a forall a. Ord a => a -> a -> Bool
< Int
2 Bool -> Bool -> Bool
|| Int
b forall a. Ord a => a -> a -> Bool
< Int
2 = (forall (v :: * -> *). Vector v Double => Int -> v Double -> Double
centralMoment Int
a v Double
xs , forall (v :: * -> *). Vector v Double => Int -> v Double -> Double
centralMoment Int
b v Double
xs)
    | Bool
otherwise      = V -> (Double, Double)
fini forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (v :: * -> *) b a.
Vector v b =>
(a -> b -> a) -> a -> v b -> a
G.foldl' V -> Double -> V
go (Double -> Double -> V
V Double
0 Double
0) forall a b. (a -> b) -> a -> b
$ v Double
xs
  where go :: V -> Double -> V
go (V Double
i Double
j) Double
x = Double -> Double -> V
V (Double
i forall a. Num a => a -> a -> a
+ Double
dDouble -> Int -> Double
^Int
a) (Double
j forall a. Num a => a -> a -> a
+ Double
dDouble -> Int -> Double
^Int
b)
            where d :: Double
d  = Double
x forall a. Num a => a -> a -> a
- Double
m
        fini :: V -> (Double, Double)
fini (V Double
i Double
j) = (Double
i forall a. Fractional a => a -> a -> a
/ Double
n , Double
j forall a. Fractional a => a -> a -> a
/ Double
n)
        m :: Double
m            = forall (v :: * -> *). Vector v Double => v Double -> Double
mean v Double
xs
        n :: Double
n            = forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v Double
xs)
{-# SPECIALIZE
    centralMoments :: Int -> Int -> U.Vector Double -> (Double, Double) #-}
{-# SPECIALIZE
    centralMoments :: Int -> Int -> V.Vector Double -> (Double, Double) #-}

-- | 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 $ U.to [1,100,101,102,103]
-- > ==> -1.497681449918257
--
-- A sample with positive skew is said to be /right-skewed/.
--
-- > skewness $ U.to [1,2,3,4,100]
-- > ==> 1.4975367033335198
--
-- A sample's skewness is not defined if its 'variance' is zero.
--
-- This function performs two passes over the sample, so is not subject
-- to stream fusion.
--
-- For samples containing many values very close to the mean, this
-- function is subject to inaccuracy due to catastrophic cancellation.
skewness :: (G.Vector v Double) => v Double -> Double
skewness :: forall (v :: * -> *). Vector v Double => v Double -> Double
skewness v Double
xs = Double
c3 forall a. Num a => a -> a -> a
* Double
c2 forall a. Floating a => a -> a -> a
** (-Double
1.5)
    where (Double
c3 , Double
c2) = forall (v :: * -> *).
Vector v Double =>
Int -> Int -> v Double -> (Double, Double)
centralMoments Int
3 Int
2 v Double
xs
{-# SPECIALIZE skewness :: U.Vector Double -> Double #-}
{-# SPECIALIZE skewness :: V.Vector Double -> Double #-}

-- | Compute the excess kurtosis of a sample.  This is a measure of
-- the \"peakedness\" of its distribution.  A high kurtosis indicates
-- that more of the sample's variance is due to infrequent severe
-- deviations, rather than more frequent modest deviations.
--
-- A sample's excess kurtosis is not defined if its 'variance' is
-- zero.
--
-- This function performs two passes over the sample, so is not subject
-- to stream fusion.
--
-- For samples containing many values very close to the mean, this
-- function is subject to inaccuracy due to catastrophic cancellation.
kurtosis :: (G.Vector v Double) => v Double -> Double
kurtosis :: forall (v :: * -> *). Vector v Double => v Double -> Double
kurtosis v Double
xs = Double
c4 forall a. Fractional a => a -> a -> a
/ (Double
c2 forall a. Num a => a -> a -> a
* Double
c2) forall a. Num a => a -> a -> a
- Double
3
    where (Double
c4 , Double
c2) = forall (v :: * -> *).
Vector v Double =>
Int -> Int -> v Double -> (Double, Double)
centralMoments Int
4 Int
2 v Double
xs
{-# SPECIALIZE kurtosis :: U.Vector Double -> Double #-}
{-# SPECIALIZE kurtosis :: V.Vector Double -> Double #-}

-- $variance
--
-- The variance — and hence the standard deviation — of a
-- sample of fewer than two elements are both defined to be zero.

-- $robust
--
-- These functions use the compensated summation algorithm of Chan et
-- al. for numerical robustness, but require two passes over the
-- sample data as a result.
--
-- Because of the need for two passes, these functions are /not/
-- subject to stream fusion.

data V = V {-# UNPACK #-} !Double {-# UNPACK #-} !Double

-- | Maximum likelihood estimate of a sample's variance.  Also known
-- as the population variance, where the denominator is /n/.
variance :: (G.Vector v Double) => v Double -> Double
variance :: forall (v :: * -> *). Vector v Double => v Double -> Double
variance v Double
samp
    | Int
n forall a. Ord a => a -> a -> Bool
> Int
1     = forall (v :: * -> *).
Vector v Double =>
Double -> v Double -> Double
robustSumVar (forall (v :: * -> *). Vector v Double => v Double -> Double
mean v Double
samp) v Double
samp forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
    | Bool
otherwise = Double
0
    where
      n :: Int
n = forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v Double
samp
{-# SPECIALIZE variance :: U.Vector Double -> Double #-}
{-# SPECIALIZE variance :: V.Vector Double -> Double #-}


-- | Unbiased estimate of a sample's variance.  Also known as the
-- sample variance, where the denominator is /n/-1.
varianceUnbiased :: (G.Vector v Double) => v Double -> Double
varianceUnbiased :: forall (v :: * -> *). Vector v Double => v Double -> Double
varianceUnbiased v Double
samp
    | Int
n forall a. Ord a => a -> a -> Bool
> Int
1     = forall (v :: * -> *).
Vector v Double =>
Double -> v Double -> Double
robustSumVar (forall (v :: * -> *). Vector v Double => v Double -> Double
mean v Double
samp) v Double
samp forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
nforall a. Num a => a -> a -> a
-Int
1)
    | Bool
otherwise = Double
0
    where
      n :: Int
n = forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v Double
samp
{-# SPECIALIZE varianceUnbiased :: U.Vector Double -> Double #-}
{-# SPECIALIZE varianceUnbiased :: V.Vector Double -> Double #-}

-- | Calculate mean and maximum likelihood estimate of variance. This
-- function should be used if both mean and variance are required
-- since it will calculate mean only once.
meanVariance ::  (G.Vector v Double) => v Double -> (Double,Double)
meanVariance :: forall (v :: * -> *).
Vector v Double =>
v Double -> (Double, Double)
meanVariance v Double
samp
  | Int
n forall a. Ord a => a -> a -> Bool
> Int
1     = (Double
m, forall (v :: * -> *).
Vector v Double =>
Double -> v Double -> Double
robustSumVar Double
m v Double
samp forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)
  | Bool
otherwise = (Double
m, Double
0)
    where
      n :: Int
n = forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v Double
samp
      m :: Double
m = forall (v :: * -> *). Vector v Double => v Double -> Double
mean v Double
samp
{-# SPECIALIZE meanVariance :: U.Vector Double -> (Double,Double) #-}
{-# SPECIALIZE meanVariance :: V.Vector Double -> (Double,Double) #-}

-- | Calculate mean and unbiased estimate of variance. This
-- function should be used if both mean and variance are required
-- since it will calculate mean only once.
meanVarianceUnb :: (G.Vector v Double) => v Double -> (Double,Double)
meanVarianceUnb :: forall (v :: * -> *).
Vector v Double =>
v Double -> (Double, Double)
meanVarianceUnb v Double
samp
  | Int
n forall a. Ord a => a -> a -> Bool
> Int
1     = (Double
m, forall (v :: * -> *).
Vector v Double =>
Double -> v Double -> Double
robustSumVar Double
m v Double
samp forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
nforall a. Num a => a -> a -> a
-Int
1))
  | Bool
otherwise = (Double
m, Double
0)
    where
      n :: Int
n = forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v Double
samp
      m :: Double
m = forall (v :: * -> *). Vector v Double => v Double -> Double
mean v Double
samp
{-# SPECIALIZE meanVarianceUnb :: U.Vector Double -> (Double,Double) #-}
{-# SPECIALIZE meanVarianceUnb :: V.Vector Double -> (Double,Double) #-}

-- | Standard deviation.  This is simply the square root of the
-- unbiased estimate of the variance.
stdDev :: (G.Vector v Double) => v Double -> Double
stdDev :: forall (v :: * -> *). Vector v Double => v Double -> Double
stdDev = forall a. Floating a => a -> a
sqrt forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (v :: * -> *). Vector v Double => v Double -> Double
varianceUnbiased
{-# SPECIALIZE stdDev :: U.Vector Double -> Double #-}
{-# SPECIALIZE stdDev :: V.Vector Double -> Double #-}

-- | Standard error of the mean. This is the standard deviation
-- divided by the square root of the sample size.
stdErrMean :: (G.Vector v Double) => v Double -> Double
stdErrMean :: forall (v :: * -> *). Vector v Double => v Double -> Double
stdErrMean v Double
samp = forall (v :: * -> *). Vector v Double => v Double -> Double
stdDev v Double
samp forall a. Fractional a => a -> a -> a
/ (forall a. Floating a => a -> a
sqrt forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (Integral a, Num b) => a -> b
fromIntegral forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (v :: * -> *) a. Vector v a => v a -> Int
G.length) v Double
samp
{-# SPECIALIZE stdErrMean :: U.Vector Double -> Double #-}
{-# SPECIALIZE stdErrMean :: V.Vector Double -> Double #-}

robustSumVarWeighted :: (G.Vector v (Double,Double)) => v (Double,Double) -> V
robustSumVarWeighted :: forall (v :: * -> *).
Vector v (Double, Double) =>
v (Double, Double) -> V
robustSumVarWeighted v (Double, Double)
samp = forall (v :: * -> *) b a.
Vector v b =>
(a -> b -> a) -> a -> v b -> a
G.foldl' V -> (Double, Double) -> V
go (Double -> Double -> V
V Double
0 Double
0) v (Double, Double)
samp
    where
      go :: V -> (Double, Double) -> V
go (V Double
s Double
w) (Double
x,Double
xw) = Double -> Double -> V
V (Double
s forall a. Num a => a -> a -> a
+ Double
xwforall a. Num a => a -> a -> a
*Double
dforall a. Num a => a -> a -> a
*Double
d) (Double
w forall a. Num a => a -> a -> a
+ Double
xw)
          where d :: Double
d = Double
x forall a. Num a => a -> a -> a
- Double
m
      m :: Double
m = forall (v :: * -> *).
Vector v (Double, Double) =>
v (Double, Double) -> Double
meanWeighted v (Double, Double)
samp
{-# INLINE robustSumVarWeighted #-}

-- | Weighted variance. This is biased estimation.
varianceWeighted :: (G.Vector v (Double,Double)) => v (Double,Double) -> Double
varianceWeighted :: forall (v :: * -> *).
Vector v (Double, Double) =>
v (Double, Double) -> Double
varianceWeighted v (Double, Double)
samp
    | forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v (Double, Double)
samp forall a. Ord a => a -> a -> Bool
> Int
1 = V -> Double
fini forall a b. (a -> b) -> a -> b
$ forall (v :: * -> *).
Vector v (Double, Double) =>
v (Double, Double) -> V
robustSumVarWeighted v (Double, Double)
samp
    | Bool
otherwise         = Double
0
    where
      fini :: V -> Double
fini (V Double
s Double
w) = Double
s forall a. Fractional a => a -> a -> a
/ Double
w
{-# SPECIALIZE varianceWeighted :: U.Vector (Double,Double) -> Double #-}
{-# SPECIALIZE varianceWeighted :: V.Vector (Double,Double) -> Double #-}

-- $cancellation
--
-- The functions prefixed with the name @fast@ below perform a single
-- pass over the sample data using Knuth's algorithm. They usually
-- work well, but see below for caveats. These functions are subject
-- to array fusion.
--
-- /Note/: in cases where most sample data is close to the sample's
-- mean, Knuth's algorithm gives inaccurate results due to
-- catastrophic cancellation.

fastVar :: (G.Vector v Double) => v Double -> T1
fastVar :: forall (v :: * -> *). Vector v Double => v Double -> T1
fastVar = forall (v :: * -> *) b a.
Vector v b =>
(a -> b -> a) -> a -> v b -> a
G.foldl' T1 -> Double -> T1
go (Int -> Double -> Double -> T1
T1 Int
0 Double
0 Double
0)
  where
    go :: T1 -> Double -> T1
go (T1 Int
n Double
m Double
s) Double
x = Int -> Double -> Double -> T1
T1 Int
n' Double
m' Double
s'
      where n' :: Int
n' = Int
n forall a. Num a => a -> a -> a
+ Int
1
            m' :: Double
m' = Double
m forall a. Num a => a -> a -> a
+ Double
d forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n'
            s' :: Double
s' = Double
s forall a. Num a => a -> a -> a
+ Double
d forall a. Num a => a -> a -> a
* (Double
x forall a. Num a => a -> a -> a
- Double
m')
            d :: Double
d  = Double
x forall a. Num a => a -> a -> a
- Double
m

-- | Maximum likelihood estimate of a sample's variance.
fastVariance :: (G.Vector v Double) => v Double -> Double
fastVariance :: forall (v :: * -> *). Vector v Double => v Double -> Double
fastVariance = T1 -> Double
fini forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (v :: * -> *). Vector v Double => v Double -> T1
fastVar
  where fini :: T1 -> Double
fini (T1 Int
n Double
_m Double
s)
          | Int
n forall a. Ord a => a -> a -> Bool
> Int
1     = Double
s forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
          | Bool
otherwise = Double
0
{-# INLINE fastVariance #-}

-- | Unbiased estimate of a sample's variance.
fastVarianceUnbiased :: (G.Vector v Double) => v Double -> Double
fastVarianceUnbiased :: forall (v :: * -> *). Vector v Double => v Double -> Double
fastVarianceUnbiased = T1 -> Double
fini forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (v :: * -> *). Vector v Double => v Double -> T1
fastVar
  where fini :: T1 -> Double
fini (T1 Int
n Double
_m Double
s)
          | Int
n forall a. Ord a => a -> a -> Bool
> Int
1     = Double
s forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
n forall a. Num a => a -> a -> a
- Int
1)
          | Bool
otherwise = Double
0
{-# INLINE fastVarianceUnbiased #-}

-- | Standard deviation.  This is simply the square root of the
-- maximum likelihood estimate of the variance.
fastStdDev :: (G.Vector v Double) => v Double -> Double
fastStdDev :: forall (v :: * -> *). Vector v Double => v Double -> Double
fastStdDev = forall a. Floating a => a -> a
sqrt forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (v :: * -> *). Vector v Double => v Double -> Double
fastVariance
{-# INLINE fastStdDev #-}

-- | Covariance of sample of pairs. For empty sample it's set to
--   zero
covariance :: (G.Vector v (Double,Double), G.Vector v Double)
           => v (Double,Double)
           -> Double
covariance :: forall (v :: * -> *).
(Vector v (Double, Double), Vector v Double) =>
v (Double, Double) -> Double
covariance v (Double, Double)
xy
  | Int
n forall a. Eq a => a -> a -> Bool
== Int
0    = Double
0
  | Bool
otherwise = forall (v :: * -> *). Vector v Double => v Double -> Double
mean forall a b. (a -> b) -> a -> b
$ forall (v :: * -> *) a b c.
(Vector v a, Vector v b, Vector v c) =>
(a -> b -> c) -> v a -> v b -> v c
G.zipWith forall a. Num a => a -> a -> a
(*)
                         (forall (v :: * -> *) a b.
(Vector v a, Vector v b) =>
(a -> b) -> v a -> v b
G.map (\Double
x -> Double
x forall a. Num a => a -> a -> a
- Double
muX) v Double
xs)
                         (forall (v :: * -> *) a b.
(Vector v a, Vector v b) =>
(a -> b) -> v a -> v b
G.map (\Double
y -> Double
y forall a. Num a => a -> a -> a
- Double
muY) v Double
ys)
  where
    n :: Int
n       = forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v (Double, Double)
xy
    (v Double
xs,v Double
ys) = forall (v :: * -> *) a b.
(Vector v a, Vector v b, Vector v (a, b)) =>
v (a, b) -> (v a, v b)
G.unzip v (Double, Double)
xy
    muX :: Double
muX     = forall (v :: * -> *). Vector v Double => v Double -> Double
mean v Double
xs
    muY :: Double
muY     = forall (v :: * -> *). Vector v Double => v Double -> Double
mean v Double
ys
{-# SPECIALIZE covariance :: U.Vector (Double,Double) -> Double #-}
{-# SPECIALIZE covariance :: V.Vector (Double,Double) -> Double #-}

-- | Correlation coefficient for sample of pairs. Also known as
--   Pearson's correlation. For empty sample it's set to zero.
correlation :: (G.Vector v (Double,Double), G.Vector v Double)
           => v (Double,Double)
           -> Double
correlation :: forall (v :: * -> *).
(Vector v (Double, Double), Vector v Double) =>
v (Double, Double) -> Double
correlation v (Double, Double)
xy
  | Int
n forall a. Eq a => a -> a -> Bool
== Int
0    = Double
0
  | Bool
otherwise = Double
cov forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt (Double
varX forall a. Num a => a -> a -> a
* Double
varY)
  where
    n :: Int
n       = forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v (Double, Double)
xy
    (v Double
xs,v Double
ys) = forall (v :: * -> *) a b.
(Vector v a, Vector v b, Vector v (a, b)) =>
v (a, b) -> (v a, v b)
G.unzip v (Double, Double)
xy
    (Double
muX,Double
varX) = forall (v :: * -> *).
Vector v Double =>
v Double -> (Double, Double)
meanVariance v Double
xs
    (Double
muY,Double
varY) = forall (v :: * -> *).
Vector v Double =>
v Double -> (Double, Double)
meanVariance v Double
ys
    cov :: Double
cov = forall (v :: * -> *). Vector v Double => v Double -> Double
mean forall a b. (a -> b) -> a -> b
$ forall (v :: * -> *) a b c.
(Vector v a, Vector v b, Vector v c) =>
(a -> b -> c) -> v a -> v b -> v c
G.zipWith forall a. Num a => a -> a -> a
(*)
            (forall (v :: * -> *) a b.
(Vector v a, Vector v b) =>
(a -> b) -> v a -> v b
G.map (\Double
x -> Double
x forall a. Num a => a -> a -> a
- Double
muX) v Double
xs)
            (forall (v :: * -> *) a b.
(Vector v a, Vector v b) =>
(a -> b) -> v a -> v b
G.map (\Double
y -> Double
y forall a. Num a => a -> a -> a
- Double
muY) v Double
ys)
{-# SPECIALIZE correlation :: U.Vector (Double,Double) -> Double #-}
{-# SPECIALIZE correlation :: V.Vector (Double,Double) -> Double #-}


-- | Pair two samples. It's like 'G.zip' but requires that both
--   samples have equal size.
pair :: (G.Vector v a, G.Vector v b, G.Vector v (a,b)) => v a -> v b -> v (a,b)
pair :: forall (v :: * -> *) a b.
(Vector v a, Vector v b, Vector v (a, b)) =>
v a -> v b -> v (a, b)
pair v a
va v b
vb
  | forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v a
va forall a. Eq a => a -> a -> Bool
== forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v b
vb = forall (v :: * -> *) a b.
(Vector v a, Vector v b, Vector v (a, b)) =>
v a -> v b -> v (a, b)
G.zip v a
va v b
vb
  | Bool
otherwise = forall a. HasCallStack => [Char] -> a
error [Char]
"Statistics.Sample.pair: vector must have same length"
{-# INLINE pair #-}

------------------------------------------------------------------------
-- Helper code. Monomorphic unpacked accumulators.

-- (^) operator from Prelude is just slow.
(^) :: Double -> Int -> Double
Double
x ^ :: Double -> Int -> Double
^ Int
1 = Double
x
Double
x ^ Int
n = Double
x forall a. Num a => a -> a -> a
* (Double
x Double -> Int -> Double
^ (Int
nforall a. Num a => a -> a -> a
-Int
1))
{-# INLINE (^) #-}

-- don't support polymorphism, as we can't get unboxed returns if we use it.
data T = T {-# UNPACK #-}!Double {-# UNPACK #-}!Int

data T1 = T1 {-# UNPACK #-}!Int {-# UNPACK #-}!Double {-# UNPACK #-}!Double

{-

Consider this core:

with data T a = T !a !Int

$wfold :: Double#
               -> Int#
               -> Int#
               -> (# Double, Int# #)

and without,

$wfold :: Double#
               -> Int#
               -> Int#
               -> (# Double#, Int# #)

yielding to boxed returns and heap checks.

-}

-- $references
--
-- * Chan, T. F.; Golub, G.H.; LeVeque, R.J. (1979) Updating formulae
--   and a pairwise algorithm for computing sample
--   variances. Technical Report STAN-CS-79-773, Department of
--   Computer Science, Stanford
--   University. <ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/79/773/CS-TR-79-773.pdf>
--
-- * Knuth, D.E. (1998) The art of computer programming, volume 2:
--   seminumerical algorithms, 3rd ed., p. 232.
--
-- * Welford, B.P. (1962) Note on a method for calculating corrected
--   sums of squares and products. /Technometrics/
--   4(3):419&#8211;420. <http://www.jstor.org/stable/1266577>
--
-- * West, D.H.D. (1979) Updating mean and variance estimates: an
--   improved method. /Communications of the ACM/
--   22(9):532&#8211;535. <http://doi.acm.org/10.1145/359146.359153>