-- |
-- Module      :  Numeric.Stats
-- Copyright   :  (c) OleksandrZhabenko 2020-2023
-- License     :  MIT
-- Stability   :  Experimental
-- Maintainer  :  oleksandr.zhabenko@yahoo.com
--
-- A very basic descriptive statistics. Functions use a tail recursion approach to compute the values and are strict by an accumulator.

{-# LANGUAGE BangPatterns, MagicHash, NoImplicitPrelude #-}

module Numeric.Stats where

import GHC.Exts
import GHC.Prim
import GHC.Base
import GHC.Real
import GHC.Float
import GHC.Num

-- | Uses GHC unlifted types from @ghc-prim@ package.
-- Similar code is here: Don Stewart.  https://donsbot.wordpress.com/2008/05/06/write-haskell-as-fast-as-c-exploiting-strictness-laziness-and-recursion/
-- And here: http://fixpt.de/blog/2017-12-04-strictness-analysis-part-1.html
-- And here: Michael Snoyman. https://www.fpcomplete.com/blog/2017/09/all-about-strictness/
--
mean2F :: [Float] -> Float# -> Int# -> Float
mean2F :: [Float] -> Float# -> Int# -> Float
mean2F ((F# !Float#
x):[Float]
xs) !Float#
s1 !Int#
l1 = [Float] -> Float# -> Int# -> Float
mean2F [Float]
xs (Float# -> Float# -> Float#
plusFloat# Float#
s1 Float#
x) (Int#
l1 Int# -> Int# -> Int#
+# Int#
1#)
mean2F [Float]
_ !Float#
s1 !Int#
l1 =
 case Int# -> Int
I# Int#
l1 of
  Int
0 -> forall a. HasCallStack => [Char] -> a
error [Char]
"Numeric.Stats.mean2F: Not defined for the third zero argument. "
  Int
_  -> Float# -> Float
F# Float#
m
   where !m :: Float#
m = Float# -> Float# -> Float#
divideFloat# Float#
s1 (Int# -> Float#
int2Float# Int#
l1)

-- | Uses GHC unlifted types from @ghc-prim@ package.
-- Similar code is here: Don Stewart.  https://donsbot.wordpress.com/2008/05/06/write-haskell-as-fast-as-c-exploiting-strictness-laziness-and-recursion/
-- And here: http://fixpt.de/blog/2017-12-04-strictness-analysis-part-1.html
-- And here: Michael Snoyman. https://www.fpcomplete.com/blog/2017/09/all-about-strictness/
--
mean2D :: [Double] -> Double# -> Int# -> Double
mean2D :: [Double] -> Double# -> Int# -> Double
mean2D ((D# !Double#
x):[Double]
xs) !Double#
s1 !Int#
l1 = [Double] -> Double# -> Int# -> Double
mean2D [Double]
xs (Double#
s1 Double# -> Double# -> Double#
+## Double#
x) (Int#
l1 Int# -> Int# -> Int#
+# Int#
1#)
mean2D [Double]
_ !Double#
s1 !Int#
l1 =
 case Int# -> Int
I# Int#
l1 of
  Int
0 -> forall a. HasCallStack => [Char] -> a
error [Char]
"Numeric.Stats.mean2D: Not defined for the third zero argument. "
  Int
_  -> Double# -> Double
D# Double#
m
   where !m :: Double#
m = Double#
s1 Double# -> Double# -> Double#
/## Int# -> Double#
int2Double# Int#
l1

-- | One-pass and tail-recursive realization for the pair of the mean and dispersion. Is vulnerable to the floating-point cancellation errors.
-- Similar code is here: Don Stewart.  https://donsbot.wordpress.com/2008/05/06/write-haskell-as-fast-as-c-exploiting-strictness-laziness-and-recursion/
-- And here: http://fixpt.de/blog/2017-12-04-strictness-analysis-part-1.html
-- And here: Michael Snoyman. https://www.fpcomplete.com/blog/2017/09/all-about-strictness/
-- When using the needed, please, refer better to their variants.
--
-- Among the 'meanWithDispersion', 'meanWithDisprsionF' and 'meanWithDispersionD' better to use the last one.
meanWithDispersionP :: (RealFrac a, Floating a) => [a] -> a -> a -> a -> a -> a -> (a,a)
meanWithDispersionP :: forall a.
(RealFrac a, Floating a) =>
[a] -> a -> a -> a -> a -> a -> (a, a)
meanWithDispersionP (!a
x:[a]
xs) !a
s1 !a
s2 !a
l1 a
m1 a
d = forall a.
(RealFrac a, Floating a) =>
[a] -> a -> a -> a -> a -> a -> (a, a)
meanWithDispersionP [a]
xs (a
s1 forall a. Num a => a -> a -> a
+ a
x) (a
s2 forall a. Num a => a -> a -> a
+ a
xforall a. Num a => a -> a -> a
*a
x) (a
l1 forall a. Num a => a -> a -> a
+ a
1) (forall {a}. Fractional a => a -> a -> a -> a
m0 a
s1 a
l1 a
x) (forall {a}. Fractional a => a -> a -> a -> a
m0 a
s2 a
l1 (a
xforall a. Num a => a -> a -> a
*a
x) forall a. Num a => a -> a -> a
- (forall {a}. Fractional a => a -> a -> a -> a
m0 a
s1 a
l1 a
x)forall a. Floating a => a -> a -> a
**a
2)
  where m0 :: a -> a -> a -> a
m0 !a
s3 !a
l2 !a
x = (a
s3 forall a. Num a => a -> a -> a
+ a
x) forall a. Fractional a => a -> a -> a
/ (a
l2 forall a. Num a => a -> a -> a
+ a
1)
meanWithDispersionP [a]
_ a
_ a
_ a
_ !a
m !a
d = (a
m,a
d)

-- | Among the 'meanWithDispersion', 'meanWithDisprsionF' and 'meanWithDispersionD' better to use the last one.
meanWithDispersionFP :: [Float] -> Float# -> Float# -> Int# -> (Float,Float)
meanWithDispersionFP :: [Float] -> Float# -> Float# -> Int# -> (Float, Float)
meanWithDispersionFP ((F# !Float#
x):[Float]
xs) !Float#
s1 !Float#
s2 !Int#
l1 = [Float] -> Float# -> Float# -> Int# -> (Float, Float)
meanWithDispersionFP [Float]
xs (Float# -> Float# -> Float#
plusFloat# Float#
s1 Float#
x) (Float# -> Float# -> Float#
plusFloat# Float#
s2 (Float# -> Float# -> Float#
timesFloat# Float#
x Float#
x)) (Int#
l1 Int# -> Int# -> Int#
+# Int#
1#)
meanWithDispersionFP [] !Float#
s1 !Float#
s2 !Int#
l1 = (Float# -> Float
F# Float#
m, Float# -> Float
F# (Float# -> Float# -> Float#
minusFloat# (Float# -> Float# -> Float#
divideFloat# Float#
s2 (Int# -> Float#
int2Float# Int#
l1)) (Float# -> Float# -> Float#
timesFloat# Float#
m Float#
m)))
  where !m :: Float#
m = Float# -> Float# -> Float#
divideFloat# Float#
s1 (Int# -> Float#
int2Float# Int#
l1)

-- | Among the 'meanWithDispersion', 'meanWithDisprsionF' and 'meanWithDispersionD' better to use the last one.
meanWithDispersionDP :: [Double] -> Double# -> Double# -> Int# -> (Double,Double)
meanWithDispersionDP :: [Double] -> Double# -> Double# -> Int# -> (Double, Double)
meanWithDispersionDP ((D# !Double#
x):[Double]
xs) !Double#
s1 !Double#
s2 !Int#
l1 = [Double] -> Double# -> Double# -> Int# -> (Double, Double)
meanWithDispersionDP [Double]
xs (Double#
s1 Double# -> Double# -> Double#
+## Double#
x) (Double#
s2 Double# -> Double# -> Double#
+## (Double#
x Double# -> Double# -> Double#
*## Double#
x)) (Int#
l1 Int# -> Int# -> Int#
+# Int#
1#)
meanWithDispersionDP [] !Double#
s1 !Double#
s2 !Int#
l1 = (Double# -> Double
D# Double#
m, Double# -> Double
D# ((Double#
s2 Double# -> Double# -> Double#
/## Int# -> Double#
int2Double# Int#
l1) Double# -> Double# -> Double#
-## (Double#
m Double# -> Double# -> Double#
*## Double#
m)))
  where !m :: Double#
m = Double#
s1 Double# -> Double# -> Double#
/## Int# -> Double#
int2Double# Int#
l1

-- | Uses 'mean2F' inside.
meanF :: [Float] -> Float
meanF :: [Float] -> Float
meanF [Float]
xs = [Float] -> Float# -> Int# -> Float
mean2F [Float]
xs Float#
0.0# Int#
0#
{-# INLINE meanF #-}

meanD :: [Double] -> Double
meanD :: [Double] -> Double
meanD [Double]
xs = [Double] -> Double# -> Int# -> Double
mean2D [Double]
xs Double#
0.0## Int#
0#

-- | Among the 'meanWithDispP', 'meanWithDispF2P' and 'meanWithDispD2P' better to use the last one.
meanWithDispP :: (RealFrac a, Floating a) => [a] -> (a,a)
meanWithDispP :: forall a. (RealFrac a, Floating a) => [a] -> (a, a)
meanWithDispP xs :: [a]
xs@(a
_:[a]
_) = forall a.
(RealFrac a, Floating a) =>
[a] -> a -> a -> a -> a -> a -> (a, a)
meanWithDispersionP [a]
xs a
0.0 a
0.0 a
0.0 a
0.0 a
0.0
meanWithDispP [a]
_ = forall a. HasCallStack => [Char] -> a
error [Char]
"Numeric.Stats.meanWithDispP: Not defined for the empty list. "
{-# RULES "realfrac/float" meanWithDispP = meanWithDispF2P #-}
{-# RULES "realfrac/double" meanWithDispP = meanWithDispD2P #-}
{-# INLINE[2] meanWithDispP #-}

-- | Among the 'meanWithDispP', 'meanWithDispF2P' and 'meanWithDispD2P' better to use the last one.
meanWithDispF2P :: [Float] -> (Float,Float)
meanWithDispF2P :: [Float] -> (Float, Float)
meanWithDispF2P xs :: [Float]
xs@(Float
_:[Float]
_) = [Float] -> Float# -> Float# -> Int# -> (Float, Float)
meanWithDispersionFP [Float]
xs Float#
0.0# Float#
0.0# Int#
0#
meanWithDispF2P [Float]
_ = forall a. HasCallStack => [Char] -> a
error [Char]
"Numeric.Stats.meanWithDispF2P: Not defined for the empty list. "
{-# INLINE meanWithDispF2P #-}

-- | Among the 'meanWithDispP', 'meanWithDispF2P' and 'meanWithDispD2P' better to use the last one.
meanWithDispD2P :: [Double] -> (Double,Double)
meanWithDispD2P :: [Double] -> (Double, Double)
meanWithDispD2P xs :: [Double]
xs@(Double
x:[Double]
_) = [Double] -> Double# -> Double# -> Int# -> (Double, Double)
meanWithDispersionDP [Double]
xs Double#
0.0## Double#
0.0## Int#
0#
meanWithDispD2P [Double]
_ = forall a. HasCallStack => [Char] -> a
error [Char]
"Numeric.Stats.meanWithDispD2P: Not defined for the empty list. "
{-# INLINE meanWithDispD2P #-}

--------------------------------------------------

-- Inspired by: https://www.khanacademy.org/math/ap-statistics/summarizing-quantitative-data-ap/more-standard-deviation/v/simulation-showing-bias-in-sample-variance
-- and:
-- https://www.khanacademy.org/math/ap-statistics/summarizing-quantitative-data-ap/more-standard-deviation/v/simulation-providing-evidence-that-n-1-gives-us-unbiased-estimate



-- | One-pass and tail-recursive realization for the pair of the mean and dispersion. Is vulnerable to the floating-point cancellation errors.
-- Similar code is here: Don Stewart.  https://donsbot.wordpress.com/2008/05/06/write-haskell-as-fast-as-c-exploiting-strictness-laziness-and-recursion/
-- And here: http://fixpt.de/blog/2017-12-04-strictness-analysis-part-1.html
-- And here: Michael Snoyman. https://www.fpcomplete.com/blog/2017/09/all-about-strictness/
-- When using the needed, please, refer better to their variants.
-- Among the 'meanWithDisp', 'meanWithDispF2' and 'meanWithDispD2' better to use the last one.
meanWithDisp :: (RealFrac a, Floating a) => [a] -> (a,a)
meanWithDisp :: forall a. (RealFrac a, Floating a) => [a] -> (a, a)
meanWithDisp xs :: [a]
xs@(a
_:a
_:[a]
_) = forall {b}. Floating b => [b] -> b -> b -> b -> (b, b)
mdl [a]
xs a
0.0 a
0.0 a
0.0
  where mdl :: [b] -> b -> b -> b -> (b, b)
mdl (!b
x:[b]
ys) !b
s1 !b
s2 !b
l1 = [b] -> b -> b -> b -> (b, b)
mdl [b]
ys (b
s1 forall a. Num a => a -> a -> a
+ b
x) (b
s2 forall a. Num a => a -> a -> a
+ b
xforall a. Num a => a -> a -> a
*b
x) (b
l1 forall a. Num a => a -> a -> a
+ b
1)
        mdl [b]
_ !b
s1 !b
s2 !b
l = (b
mm, (b
s2 forall a. Num a => a -> a -> a
- b
mmforall a. Floating a => a -> a -> a
**b
2 forall a. Num a => a -> a -> a
* b
l) forall a. Fractional a => a -> a -> a
/ (b
l forall a. Num a => a -> a -> a
- b
1))
          where !mm :: b
mm = b
s1 forall a. Fractional a => a -> a -> a
/ b
l
meanWithDisp [a]
_ = forall a. HasCallStack => [Char] -> a
error [Char]
"Numeric.Stats.meanWithDisp: Not defined for the list with less than two elements. "
{-# RULES "realfrac/float" meanWithDisp = meanWithDispF2 #-}
{-# RULES "realfrac/double" meanWithDisp = meanWithDispD2 #-}
{-# INLINE[2] meanWithDisp #-}

-- | One-pass and tail-recursive realization for the pair of the mean and dispersion. Is vulnerable to the floating-point cancellation errors.
-- Similar code is here: Don Stewart.  https://donsbot.wordpress.com/2008/05/06/write-haskell-as-fast-as-c-exploiting-strictness-laziness-and-recursion/
-- And here: http://fixpt.de/blog/2017-12-04-strictness-analysis-part-1.html
-- And here: Michael Snoyman. https://www.fpcomplete.com/blog/2017/09/all-about-strictness/
-- When using the needed, please, refer better to their variants.
-- Among the 'meanWithDisp', 'meanWithDispF2' and 'meanWithDispD2' better to use the last one.
meanWithDispF2 :: [Float] -> (Float,Float)
meanWithDispF2 :: [Float] -> (Float, Float)
meanWithDispF2 xs :: [Float]
xs@(Float
_:Float
_:[Float]
_) = [Float] -> Float# -> Float# -> Int# -> (Float, Float)
mdlF [Float]
xs Float#
0.0# Float#
0.0# Int#
0#
  where mdlF :: [Float] -> Float# -> Float# -> Int# -> (Float, Float)
mdlF (F# !Float#
x:[Float]
ys) !Float#
s1 !Float#
s2 !Int#
l1 = [Float] -> Float# -> Float# -> Int# -> (Float, Float)
mdlF [Float]
ys (Float# -> Float# -> Float#
plusFloat# Float#
s1 Float#
x) (Float# -> Float# -> Float#
plusFloat# Float#
s2 (Float# -> Float# -> Float#
timesFloat# Float#
x Float#
x)) (Int#
l1 Int# -> Int# -> Int#
+# Int#
1#)
        mdlF [] !Float#
s1 !Float#
s2 !Int#
l1 = (Float# -> Float
F# Float#
m, Float# -> Float
F# (Float# -> Float# -> Float#
divideFloat# (Float# -> Float# -> Float#
minusFloat# Float#
s2 (Float# -> Float# -> Float#
timesFloat# (Float# -> Float# -> Float#
timesFloat# Float#
m Float#
m) (Int# -> Float#
int2Float# Int#
l1))) (Int# -> Float#
int2Float# (Int#
l1 Int# -> Int# -> Int#
-# Int#
1#))))
          where !m :: Float#
m = Float# -> Float# -> Float#
divideFloat# Float#
s1 (Int# -> Float#
int2Float# Int#
l1)
meanWithDispF2 [Float]
_ = forall a. HasCallStack => [Char] -> a
error [Char]
"Numeric.Stats.meanWithDispF2: Not defined for the list with less than two elements. "
{-# INLINE meanWithDispF2 #-}

-- | One-pass and tail-recursive realization for the pair of the mean and dispersion. Is vulnerable to the floating-point cancellation errors.
-- Similar code is here: Don Stewart.  https://donsbot.wordpress.com/2008/05/06/write-haskell-as-fast-as-c-exploiting-strictness-laziness-and-recursion/
-- And here: http://fixpt.de/blog/2017-12-04-strictness-analysis-part-1.html
-- And here: Michael Snoyman. https://www.fpcomplete.com/blog/2017/09/all-about-strictness/
-- When using the needed, please, refer better to their variants.
-- Among the 'meanWithDisp', 'meanWithDispF2' and 'meanWithDispD2' better to use the last one.
meanWithDispD2 :: [Double] -> (Double,Double)
meanWithDispD2 :: [Double] -> (Double, Double)
meanWithDispD2 xs :: [Double]
xs@(Double
_:Double
_:[Double]
_) = [Double] -> Double# -> Double# -> Int# -> (Double, Double)
mdlD [Double]
xs Double#
0.0## Double#
0.0## Int#
0#
  where mdlD :: [Double] -> Double# -> Double# -> Int# -> (Double, Double)
mdlD ((D# !Double#
x):[Double]
xs) !Double#
s1 !Double#
s2 !Int#
l1 = [Double] -> Double# -> Double# -> Int# -> (Double, Double)
mdlD [Double]
xs (Double#
s1 Double# -> Double# -> Double#
+## Double#
x) (Double#
s2 Double# -> Double# -> Double#
+## (Double#
x Double# -> Double# -> Double#
*## Double#
x)) (Int#
l1 Int# -> Int# -> Int#
+# Int#
1#)
        mdlD [] !Double#
s1 !Double#
s2 !Int#
l1 = (Double# -> Double
D# Double#
m, Double# -> Double
D# ((Double#
s2 Double# -> Double# -> Double#
-## Double#
m Double# -> Double# -> Double#
*## Double#
m Double# -> Double# -> Double#
*## Int# -> Double#
int2Double# Int#
l1) Double# -> Double# -> Double#
/## (Int# -> Double#
int2Double# (Int#
l1 Int# -> Int# -> Int#
-# Int#
1#))))
          where !m :: Double#
m = Double#
s1 Double# -> Double# -> Double#
/## Int# -> Double#
int2Double# Int#
l1
meanWithDispD2 [Double]
_ = forall a. HasCallStack => [Char] -> a
error [Char]
"Numeric.Stats.meanWithDispD2: Not defined for the list with less than two elements. "
{-# INLINE meanWithDispD2 #-}