-- SPDX-License-Identifier: MPL-2.0

{- | Finite-data checked statistics. Quantiles use Hyndman--Fan type 7; sample
variance divides by n-1; skewness and excess kurtosis use population central
moments. Weights are reliability weights normalized by their positive sum.
NaN/infinity are rejected; 'masked' selects True entries before validation.
Hypothesis tests report a statistic, a two-sided p value (unless intrinsically
upper-tailed), and degrees of freedom. Rank tests use tie-corrected normal/
chi-square asymptotics. KS uses the asymptotic continuous-null law; normality
tests require a fully specified or fitted normal as stated on each function.
-}
module HMath.Statistics.Advanced (
    arithmeticMean,
    weightedMean,
    geometricMean,
    harmonicMean,
    rms,
    trimmedMean,
    winsorizedMean,
    median,
    modes,
    quantile,
    percentile,
    interquartileRange,
    medianAbsoluteDeviation,
    dataRange,
    populationVariance,
    sampleVariance,
    populationSD,
    sampleSD,
    coefficientOfVariation,
    centralMoment,
    skewness,
    kurtosis,
    covariance,
    correlation,
    covarianceMatrix,
    spearman,
    kendall,
    partialCorrelation,
    autocorrelation,
    crossCorrelation,
    movingAverage,
    exponentialMovingAverage,
    movingVariance,
    Streaming,
    emptyStreaming,
    updateStreaming,
    combineStreaming,
    streamingResult,
    histogramND,
    entropy,
    mutualInformation,
    klDivergence,
    jsDivergence,
    zScores,
    standardError,
    confidenceInterval,
    bootstrap,
    jackknife,
    TestResult (..),
    oneSampleT,
    twoSampleT,
    pairedT,
    zTest,
    chiSquareTest,
    varianceFTest,
    anova,
    welchAnova,
    kruskalWallis,
    mannWhitney,
    wilcoxon,
    kolmogorovSmirnov,
    shapiroWilk,
    andersonDarling,
    levene,
    bartlett,
    fisherExact,
    permutationTest,
    bonferroni,
    benjaminiHochberg,
    cohensD,
    etaSquared,
    normalTestPower,
    normalSampleSize,
    tukeyOutliers,
    grubbs,
    masked,
    huberLocation,
) where

import qualified Data.Foldable as Foldable (foldl')
import Data.List (group, sort, sortOn, transpose)
import HMath.Error
import HMath.Matrix (fromListChecked, toList)
import HMath.Matrix.LU (inverse)
import qualified HMath.Probability as P
import qualified HMath.Random as R
import qualified HMath.Special as S

finite :: Double -> Bool
finite :: Double -> Bool
finite Double
x = Bool -> Bool
not (Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x)
require :: Bool -> String -> Either MathError ()
require :: Bool -> String -> Either MathError ()
require Bool
p String
s = if Bool
p then () -> Either MathError ()
forall a b. b -> Either a b
Right () else MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
s)
checked :: Double -> Either MathError Double
checked :: Double -> Either MathError Double
checked Double
x = if Double -> Bool
finite Double
x then Double -> Either MathError Double
forall a b. b -> Either a b
Right Double
x else MathError -> Either MathError Double
forall a b. a -> Either a b
Left MathError
NonFiniteResult
validate :: Int -> [Double] -> Either MathError ()
validate :: Int -> [Double] -> Either MathError ()
validate Int
n [Double]
xs =
    Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
n Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
xs) (String
"requires at least " String -> String -> String
forall a. [a] -> [a] -> [a]
++ Int -> String
forall a. Show a => a -> String
show Int
n String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
" finite observations")

-- | Mean of finite nonempty data. Dividing each observation before summation avoids preventable same-sign sum overflow.
arithmeticMean :: [Double] -> Either MathError Double
arithmeticMean :: [Double] -> Either MathError Double
arithmeticMean [Double]
xs = do Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs; Double -> Either MathError Double
checked ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs)) [Double]
xs))

-- | Reliability-weighted mean of observations and nonnegative finite weights; weights must have finite positive total.
weightedMean :: [Double] -> [Double] -> Either MathError Double
weightedMean :: [Double] -> [Double] -> Either MathError Double
weightedMean [Double]
xs [Double]
ws = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs
    Bool -> String -> Either MathError ()
require
        ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ws Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Double
w -> Double -> Bool
finite Double
w Bool -> Bool -> Bool
&& Double
w Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0) [Double]
ws Bool -> Bool -> Bool
&& Double -> Bool
finite ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
ws) Bool -> Bool -> Bool
&& [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
ws Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0)
        String
"invalid reliability weights"
    Double -> Either MathError Double
checked ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
x Double
w -> Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
w Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
ws)) [Double]
xs [Double]
ws))

-- | exp(mean(log x)) for strictly positive finite observations.
geometricMean :: [Double] -> Either MathError Double
geometricMean :: [Double] -> Either MathError Double
geometricMean [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs
    Bool -> String -> Either MathError ()
require ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) [Double]
xs) String
"geometric mean requires positive data"
    m <- [Double] -> Either MathError Double
arithmeticMean ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Floating a => a -> a
log [Double]
xs)
    checked (exp m)

-- | Reciprocal of the mean reciprocal, for strictly positive finite observations.
harmonicMean :: [Double] -> Either MathError Double
harmonicMean :: [Double] -> Either MathError Double
harmonicMean [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs
    Bool -> String -> Either MathError ()
require ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) [Double]
xs) String
"harmonic mean requires positive data"
    m <- [Double] -> Either MathError Double
arithmeticMean ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/) [Double]
xs)
    checked (1 / m)

-- | Root-mean-square with max scaling to avoid squaring overflow.
rms :: [Double] -> Either MathError Double
rms :: [Double] -> Either MathError Double
rms [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs
    let scale :: Double
scale = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs [Double]
xs)
     in if Double
scale Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0
            then Double -> Either MathError Double
forall a. a -> Either MathError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure Double
0
            else
                Double -> Either MathError Double
checked (Double
scale Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Double
x -> (Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
scale) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) [Double]
xs) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs)))

-- | Discard floor(fraction*n) samples from each sorted tail, with fraction in [0,1/2), then average.
trimmedMean :: Double -> [Double] -> Either MathError Double
trimmedMean :: Double -> [Double] -> Either MathError Double
trimmedMean Double
fraction [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
fraction Bool -> Bool -> Bool
&& Double
fraction Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 Bool -> Bool -> Bool
&& Double
fraction Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0.5) String
"trim fraction must lie in [0,0.5)"
    let k :: Int
k = Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double
fraction Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs))
     in [Double] -> Either MathError Double
arithmeticMean (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
k ([Double] -> [Double]
forall a. Ord a => [a] -> [a]
sort [Double]
xs)))

-- | Replace floor(fraction*n) values at each tail with the nearest retained value, then average.
winsorizedMean :: Double -> [Double] -> Either MathError Double
winsorizedMean :: Double -> [Double] -> Either MathError Double
winsorizedMean Double
fraction [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
fraction Bool -> Bool -> Bool
&& Double
fraction Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 Bool -> Bool -> Bool
&& Double
fraction Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0.5) String
"winsor fraction must lie in [0,0.5)"
    let ys :: [Double]
ys = [Double] -> [Double]
forall a. Ord a => [a] -> [a]
sort [Double]
xs
        k :: Int
k = Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double
fraction Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs))
        lo :: Double
lo = [Double]
ys [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
k
        hi :: Double
hi = [Double]
ys [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
     in [Double] -> Either MathError Double
arithmeticMean ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
lo (Double -> Double) -> (Double -> Double) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
hi) [Double]
xs)

-- | Hyndman--Fan type-7 quantile with linear interpolation and probability p in [0,1].
quantile :: Double -> [Double] -> Either MathError Double
quantile :: Double -> [Double] -> Either MathError Double
quantile Double
p [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
p Bool -> Bool -> Bool
&& Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 Bool -> Bool -> Bool
&& Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1) String
"quantile probability outside [0,1]"
    let ys :: [Double]
ys = [Double] -> [Double]
forall a. Ord a => [a] -> [a]
sort [Double]
xs
        t :: Double
t = Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
        i :: Int
i = Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor Double
t
        j :: Int
j = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
        w :: Double
w = Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i
     in Double -> Either MathError Double
checked ((Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
w) Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
ys [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
ys [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j))

-- | Type-7 50th percentile; averages the central pair for an even sample count.
median :: [Double] -> Either MathError Double
median :: [Double] -> Either MathError Double
median = Double -> [Double] -> Either MathError Double
quantile Double
0.5

-- | Type-7 quantile with percentage in [0,100].
percentile :: Double -> [Double] -> Either MathError Double
percentile :: Double -> [Double] -> Either MathError Double
percentile Double
p = Double -> [Double] -> Either MathError Double
quantile (Double
p Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
100)

-- | All values of maximum multiplicity, sorted ascending; exact floating equality defines ties.
modes :: [Double] -> Either MathError [Double]
modes :: [Double] -> Either MathError [Double]
modes [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs
    let gs :: [[Double]]
gs = [Double] -> [[Double]]
forall a. Eq a => [a] -> [[a]]
group ([Double] -> [Double]
forall a. Ord a => [a] -> [a]
sort [Double]
xs); m :: Int
m = [Int] -> Int
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (([Double] -> Int) -> [[Double]] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
gs) in [Double] -> Either MathError [Double]
forall a. a -> Either MathError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure [[Double]
ys [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0 | [Double]
ys <- [[Double]]
gs, [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
m]

-- | Type-7 75th percentile minus 25th percentile.
interquartileRange :: [Double] -> Either MathError Double
interquartileRange :: [Double] -> Either MathError Double
interquartileRange [Double]
xs = do a <- Double -> [Double] -> Either MathError Double
quantile Double
0.25 [Double]
xs; b <- quantile 0.75 xs; checked (b - a)

-- | Median absolute deviation from the median, without the normal consistency multiplier.
medianAbsoluteDeviation :: [Double] -> Either MathError Double
medianAbsoluteDeviation :: [Double] -> Either MathError Double
medianAbsoluteDeviation [Double]
xs = do m <- [Double] -> Either MathError Double
median [Double]
xs; median (map (\Double
x -> Double -> Double
forall a. Num a => a -> a
abs (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
m)) xs)

-- | Maximum minus minimum; arithmetic overflow is an explicit error.
dataRange :: [Double] -> Either MathError Double
dataRange :: [Double] -> Either MathError Double
dataRange [Double]
xs = do Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs; Double -> Either MathError Double
checked ([Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
xs Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
xs)

-- | Population central moment of integer order 0..32; nonfinite intermediate/final arithmetic is rejected.
centralMoment :: Int -> [Double] -> Either MathError Double
centralMoment :: Int -> [Double] -> Either MathError Double
centralMoment Int
k [Double]
xs = do
    Bool -> String -> Either MathError ()
require (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
<= Int
32) String
"moment order must be 0..32"
    m <- [Double] -> Either MathError Double
arithmeticMean [Double]
xs
    arithmeticMean (map (\Double
x -> (Double
x 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
k) xs)

-- | Second population central moment, dividing by n.
populationVariance :: [Double] -> Either MathError Double
populationVariance :: [Double] -> Either MathError Double
populationVariance = Int -> [Double] -> Either MathError Double
centralMoment Int
2

-- | Unbiased sample variance, dividing by n-1; requires at least two observations.
sampleVariance :: [Double] -> Either MathError Double
sampleVariance :: [Double] -> Either MathError Double
sampleVariance [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
2 [Double]
xs
    v <- [Double] -> Either MathError Double
populationVariance [Double]
xs
    checked (v * fromIntegral (length xs) / fromIntegral (length xs - 1))

-- | Square root of population variance.
populationSD :: [Double] -> Either MathError Double
populationSD :: [Double] -> Either MathError Double
populationSD = (Double -> Double)
-> Either MathError Double -> Either MathError Double
forall a b. (a -> b) -> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap Double -> Double
forall a. Floating a => a -> a
sqrt (Either MathError Double -> Either MathError Double)
-> ([Double] -> Either MathError Double)
-> [Double]
-> Either MathError Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Either MathError Double
populationVariance

-- | Square root of unbiased sample variance.
sampleSD :: [Double] -> Either MathError Double
sampleSD :: [Double] -> Either MathError Double
sampleSD = (Double -> Double)
-> Either MathError Double -> Either MathError Double
forall a b. (a -> b) -> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap Double -> Double
forall a. Floating a => a -> a
sqrt (Either MathError Double -> Either MathError Double)
-> ([Double] -> Either MathError Double)
-> [Double]
-> Either MathError Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Either MathError Double
sampleVariance

-- | Population standard deviation divided by absolute mean; rejects zero mean.
coefficientOfVariation :: [Double] -> Either MathError Double
coefficientOfVariation :: [Double] -> Either MathError Double
coefficientOfVariation [Double]
xs = do
    m <- [Double] -> Either MathError Double
arithmeticMean [Double]
xs
    require (m /= 0) "coefficient of variation undefined at zero mean"
    s <- populationSD xs
    checked (s / abs m)

-- | Uncorrected population third standardized central moment; requires nonconstant observations.
skewness :: [Double] -> Either MathError Double
skewness :: [Double] -> Either MathError Double
skewness [Double]
xs = do
    v <- [Double] -> Either MathError Double
populationVariance [Double]
xs
    require (v > 0) "skewness needs nonconstant data"
    m <- centralMoment 3 xs
    checked (m / v ** 1.5)

-- | Uncorrected population excess kurtosis, subtracting three; requires nonconstant observations.
kurtosis :: [Double] -> Either MathError Double
kurtosis :: [Double] -> Either MathError Double
kurtosis [Double]
xs = do
    v <- [Double] -> Either MathError Double
populationVariance [Double]
xs
    require (v > 0) "kurtosis needs nonconstant data"
    m <- centralMoment 4 xs
    checked (m / (v * v) - 3)

-- | Covariance of two matching finite data vectors; True selects the n-1 sample correction.
covariance :: Bool -> [Double] -> [Double] -> Either MathError Double
covariance :: Bool -> [Double] -> [Double] -> Either MathError Double
covariance Bool
sampleCorrection [Double]
xs [Double]
ys = do
    Int -> [Double] -> Either MathError ()
validate (if Bool
sampleCorrection then Int
2 else Int
1) [Double]
xs
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
ys
    Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys) String
"covariance length mismatch"
    a <- [Double] -> Either MathError Double
arithmeticMean [Double]
xs
    b <- arithmeticMean ys
    checked
        ( sum (zipWith (\Double
x Double
y -> (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b)) xs ys)
            / fromIntegral (length xs - if sampleCorrection then 1 else 0)
        )

-- | Pearson correlation of matching nonconstant data vectors.
correlation :: [Double] -> [Double] -> Either MathError Double
correlation :: [Double] -> [Double] -> Either MathError Double
correlation [Double]
xs [Double]
ys = do
    c <- Bool -> [Double] -> [Double] -> Either MathError Double
covariance Bool
False [Double]
xs [Double]
ys
    a <- populationSD xs
    b <- populationSD ys
    require (a > 0 && b > 0) "correlation needs nonconstant data"
    checked (max (-1) (min 1 (c / (a * b))))

-- | Covariance matrix for variables supplied as columns (each list is one variable). Boolean selects sample correction.
covarianceMatrix :: Bool -> [[Double]] -> Either MathError [[Double]]
covarianceMatrix :: Bool -> [[Double]] -> Either MathError [[Double]]
covarianceMatrix Bool
correction [[Double]]
columns = do
    Bool -> String -> Either MathError ()
require (Bool -> Bool
not ([[Double]] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [[Double]]
columns)) String
"covariance matrix needs columns"
    ([Double] -> Either MathError [Double])
-> [[Double]] -> Either MathError [[Double]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM (\[Double]
x -> ([Double] -> Either MathError Double)
-> [[Double]] -> Either MathError [Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM (Bool -> [Double] -> [Double] -> Either MathError Double
covariance Bool
correction [Double]
x) [[Double]]
columns) [[Double]]
columns

ranks :: [Double] -> [Double]
ranks :: [Double] -> [Double]
ranks [Double]
xs = ((Int, Double) -> Double) -> [(Int, Double)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Double) -> Double
forall a b. (a, b) -> b
snd (((Int, Double) -> Int) -> [(Int, Double)] -> [(Int, Double)]
forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn (Int, Double) -> Int
forall a b. (a, b) -> a
fst (([(Int, Double)] -> [(Int, Double)])
-> [[(Int, Double)]] -> [(Int, Double)]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap [(Int, Double)] -> [(Int, Double)]
forall {b} {a}. Fractional b => [(a, Double)] -> [(a, b)]
assign ([(Int, Double)] -> [[(Int, Double)]]
forall {b} {a}. Eq b => [(a, b)] -> [[(a, b)]]
groupByValue (((Int, Double) -> Double) -> [(Int, Double)] -> [(Int, Double)]
forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn (Int, Double) -> Double
forall a b. (a, b) -> b
snd ([Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [Double]
xs)))))
  where
    groupByValue :: [(a, b)] -> [[(a, b)]]
groupByValue [] = []
    groupByValue ((a, b)
a : [(a, b)]
as) = let ([(a, b)]
g, [(a, b)]
t) = ((a, b) -> Bool) -> [(a, b)] -> ([(a, b)], [(a, b)])
forall a. (a -> Bool) -> [a] -> ([a], [a])
span ((b -> b -> Bool
forall a. Eq a => a -> a -> Bool
== (a, b) -> b
forall a b. (a, b) -> b
snd (a, b)
a) (b -> Bool) -> ((a, b) -> b) -> (a, b) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a, b) -> b
forall a b. (a, b) -> b
snd) [(a, b)]
as in ((a, b)
a (a, b) -> [(a, b)] -> [(a, b)]
forall a. a -> [a] -> [a]
: [(a, b)]
g) [(a, b)] -> [[(a, b)]] -> [[(a, b)]]
forall a. a -> [a] -> [a]
: [(a, b)] -> [[(a, b)]]
groupByValue [(a, b)]
t
    assign :: [(a, Double)] -> [(a, b)]
assign [(a, Double)]
g =
        let v :: b
v =
                Int -> b
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Int] -> Int
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((Double -> Bool) -> [Double] -> [Double]
forall a. (a -> Bool) -> [a] -> [a]
filter (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
x) [Double]
xs) | (a
_, Double
x) <- [(a, Double)]
g]) b -> b -> b
forall a. Fractional a => a -> a -> a
/ Int -> b
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([(a, Double)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(a, Double)]
g)
                    b -> b -> b
forall a. Num a => a -> a -> a
+ Int -> b
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([(a, Double)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(a, Double)]
g Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) b -> b -> b
forall a. Fractional a => a -> a -> a
/ b
2
         in [(a
i, b
v) | (a
i, Double
_) <- [(a, Double)]
g]

-- | Pearson correlation of average tied ranks; matching nonconstant rank vectors required.
spearman :: [Double] -> [Double] -> Either MathError Double
spearman :: [Double] -> [Double] -> Either MathError Double
spearman [Double]
xs [Double]
ys = do Int -> [Double] -> Either MathError ()
validate Int
2 [Double]
xs; Int -> [Double] -> Either MathError ()
validate Int
2 [Double]
ys; [Double] -> [Double] -> Either MathError Double
correlation ([Double] -> [Double]
ranks [Double]
xs) ([Double] -> [Double]
ranks [Double]
ys)

-- | Kendall tau-b with tied-pair correction; quadratic reference implementation.
kendall :: [Double] -> [Double] -> Either MathError Double
kendall :: [Double] -> [Double] -> Either MathError Double
kendall [Double]
xs [Double]
ys = do
    Int -> [Double] -> Either MathError ()
validate Int
2 [Double]
xs
    Int -> [Double] -> Either MathError ()
validate Int
2 [Double]
ys
    Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys) String
"Kendall length mismatch"
    let pairs :: [(Ordering, Ordering)]
pairs =
            [ (Double -> Double -> Ordering
forall a. Ord a => a -> a -> Ordering
compare ([Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) ([Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j), Double -> Double -> Ordering
forall a. Ord a => a -> a -> Ordering
compare ([Double]
ys [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) ([Double]
ys [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j))
            | Int
i <- [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
            , Int
j <- [Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
            ]
        score :: Int
score = [Int] -> Int
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [if Ordering
a Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
== Ordering
b then Int
1 else -Int
1 | (Ordering
a, Ordering
b) <- [(Ordering, Ordering)]
pairs, Ordering
a Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
/= Ordering
EQ, Ordering
b Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
/= Ordering
EQ]
        n :: Double
n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([(Ordering, Ordering)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(Ordering, Ordering)]
pairs)
        tx :: Double
tx = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([(Ordering, Ordering)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (((Ordering, Ordering) -> Bool)
-> [(Ordering, Ordering)] -> [(Ordering, Ordering)]
forall a. (a -> Bool) -> [a] -> [a]
filter ((Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
== Ordering
EQ) (Ordering -> Bool)
-> ((Ordering, Ordering) -> Ordering)
-> (Ordering, Ordering)
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Ordering, Ordering) -> Ordering
forall a b. (a, b) -> a
fst) [(Ordering, Ordering)]
pairs))
        ty :: Double
ty = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([(Ordering, Ordering)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (((Ordering, Ordering) -> Bool)
-> [(Ordering, Ordering)] -> [(Ordering, Ordering)]
forall a. (a -> Bool) -> [a] -> [a]
filter ((Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
== Ordering
EQ) (Ordering -> Bool)
-> ((Ordering, Ordering) -> Ordering)
-> (Ordering, Ordering)
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Ordering, Ordering) -> Ordering
forall a b. (a, b) -> b
snd) [(Ordering, Ordering)]
pairs))
        den :: Double
den = Double -> Double
forall a. Floating a => a -> a
sqrt ((Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
tx) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
ty))
    Bool -> String -> Either MathError ()
require (Double
den Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"Kendall undefined for constant ranks"
    Double -> Either MathError Double
forall a. a -> Either MathError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
score :: Int) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
den)

-- | Partial correlation of variable indices i,j conditional on every other supplied column, through the inverse sample covariance.
partialCorrelation :: [[Double]] -> Int -> Int -> Either MathError Double
partialCorrelation :: [[Double]] -> Int -> Int -> Either MathError Double
partialCorrelation [[Double]]
columns Int
i Int
j = do
    let p :: Int
p = [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
columns
    Bool -> String -> Either MathError ()
require (Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
p Bool -> Bool -> Bool
&& Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
p Bool -> Bool -> Bool
&& Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
j) String
"partial correlation indices invalid"
    cov <- Bool -> [[Double]] -> Either MathError [[Double]]
covarianceMatrix Bool
True [[Double]]
columns
    mat <- fromListChecked p p (concat cov)
    inv <- inverse mat
    let a = Matrix Double -> [Double]
forall a. Matrix a -> [a]
toList Matrix Double
inv; den = Double -> Double
forall a. Floating a => a -> a
sqrt ([Double]
a [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
i) Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
a [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j))
    require (den > 0) "partial correlation invalid precision matrix"
    checked (-a !! (i * p + j) / den)

-- | Pearson correlation between overlapping lagged segments; zero<=lag<n-1.
autocorrelation :: Int -> [Double] -> Either MathError Double
autocorrelation :: Int -> [Double] -> Either MathError Double
autocorrelation Int
lag [Double]
xs = Int -> [Double] -> [Double] -> Either MathError Double
crossCorrelation Int
lag [Double]
xs [Double]
xs

-- | Pearson correlation of overlapping x_t and y_(t+lag), using their common finite length.
crossCorrelation :: Int -> [Double] -> [Double] -> Either MathError Double
crossCorrelation :: Int -> [Double] -> [Double] -> Either MathError Double
crossCorrelation Int
lag [Double]
xs [Double]
ys = do
    Bool -> String -> Either MathError ()
require (Int
lag Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
lag Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int -> Int -> Int
forall a. Ord a => a -> a -> a
min ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs) ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) String
"lag outside overlap"
    [Double] -> [Double] -> Either MathError Double
correlation
        (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take (Int -> Int -> Int
forall a. Ord a => a -> a -> a
min ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs) ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
lag) [Double]
xs)
        (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
lag (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take (Int -> Int -> Int
forall a. Ord a => a -> a -> a
min ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs) ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys)) [Double]
ys))

-- | Valid-window moving means; window must lie in [1,n], with no endpoint padding.
movingAverage :: Int -> [Double] -> Either MathError [Double]
movingAverage :: Int -> [Double] -> Either MathError [Double]
movingAverage Int
window [Double]
xs = ([Double] -> Either MathError Double)
-> Int -> [Double] -> Either MathError [Double]
rolling [Double] -> Either MathError Double
arithmeticMean Int
window [Double]
xs

-- | Valid-window population variances with no endpoint padding.
movingVariance :: Int -> [Double] -> Either MathError [Double]
movingVariance :: Int -> [Double] -> Either MathError [Double]
movingVariance Int
window [Double]
xs = ([Double] -> Either MathError Double)
-> Int -> [Double] -> Either MathError [Double]
rolling [Double] -> Either MathError Double
populationVariance Int
window [Double]
xs

rolling :: ([Double] -> Either MathError Double) -> Int -> [Double] -> Either MathError [Double]
rolling :: ([Double] -> Either MathError Double)
-> Int -> [Double] -> Either MathError [Double]
rolling [Double] -> Either MathError Double
f Int
n [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs
    Bool -> String -> Either MathError ()
require (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs) String
"rolling window outside data"
    (Int -> Either MathError Double)
-> [Int] -> Either MathError [Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM (\Int
i -> [Double] -> Either MathError Double
f (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
n (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
i [Double]
xs))) [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
n]

-- | EMA with alpha in (0,1], initialized at the first observation.
exponentialMovingAverage :: Double -> [Double] -> Either MathError [Double]
exponentialMovingAverage :: Double -> [Double] -> Either MathError [Double]
exponentialMovingAverage Double
alpha [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
alpha Bool -> Bool -> Bool
&& Double
alpha Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
alpha Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1) String
"EMA alpha in (0,1] required"
    [Double] -> Either MathError [Double]
forall a. a -> Either MathError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure ((Double -> Double -> Double) -> Double -> [Double] -> [Double]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\Double
a Double
x -> (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
alpha) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x) ([Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
xs))

-- | Opaque Welford accumulator storing an unbounded Integer count, mean and sum of squared deviations.
data Streaming = Streaming !Integer !Double !Double deriving (Streaming -> Streaming -> Bool
(Streaming -> Streaming -> Bool)
-> (Streaming -> Streaming -> Bool) -> Eq Streaming
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Streaming -> Streaming -> Bool
== :: Streaming -> Streaming -> Bool
$c/= :: Streaming -> Streaming -> Bool
/= :: Streaming -> Streaming -> Bool
Eq, Int -> Streaming -> String -> String
[Streaming] -> String -> String
Streaming -> String
(Int -> Streaming -> String -> String)
-> (Streaming -> String)
-> ([Streaming] -> String -> String)
-> Show Streaming
forall a.
(Int -> a -> String -> String)
-> (a -> String) -> ([a] -> String -> String) -> Show a
$cshowsPrec :: Int -> Streaming -> String -> String
showsPrec :: Int -> Streaming -> String -> String
$cshow :: Streaming -> String
show :: Streaming -> String
$cshowList :: [Streaming] -> String -> String
showList :: [Streaming] -> String -> String
Show)

-- | Zero-count Welford accumulator.
emptyStreaming :: Streaming
emptyStreaming :: Streaming
emptyStreaming = Integer -> Double -> Double -> Streaming
Streaming Integer
0 Double
0 Double
0

-- | Add one finite observation by Welford recurrence; nonfinite arithmetic is rejected.
updateStreaming :: Streaming -> Double -> Either MathError Streaming
updateStreaming :: Streaming -> Double -> Either MathError Streaming
updateStreaming (Streaming Integer
n Double
m Double
m2) Double
x = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double
x]
    let count :: Integer
count = Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1; d :: Double
d = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
m; mean :: Double
mean = Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Integer -> Double
forall a. Num a => Integer -> a
fromInteger Integer
count; second :: Double
second = Double
m2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
mean)
     in Double -> Either MathError Double
checked Double
mean Either MathError Double
-> Either MathError Double -> Either MathError Double
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Double -> Either MathError Double
checked Double
second Either MathError Double
-> Either MathError Streaming -> Either MathError Streaming
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Streaming -> Either MathError Streaming
forall a. a -> Either MathError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Integer -> Double -> Double -> Streaming
Streaming Integer
count Double
mean Double
second)

-- | Chan parallel merge of two Welford accumulators; empty accumulators are identities.
combineStreaming :: Streaming -> Streaming -> Either MathError Streaming
combineStreaming :: Streaming -> Streaming -> Either MathError Streaming
combineStreaming a :: Streaming
a@(Streaming Integer
n Double
m Double
v) b :: Streaming
b@(Streaming Integer
k Double
u Double
w)
    | Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 = Streaming -> Either MathError Streaming
forall a b. b -> Either a b
Right Streaming
b
    | Integer
k Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 = Streaming -> Either MathError Streaming
forall a b. b -> Either a b
Right Streaming
a
    | Bool
otherwise =
        let d :: Double
d = Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
m
            count :: Integer
count = Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
k
            mean :: Double
mean = Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Integer -> Double
forall a. Num a => Integer -> a
fromInteger Integer
k Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Integer -> Double
forall a. Num a => Integer -> a
fromInteger Integer
count
            second :: Double
second = Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Integer -> Double
forall a. Num a => Integer -> a
fromInteger Integer
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* Integer -> Double
forall a. Num a => Integer -> a
fromInteger Integer
k Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Integer -> Double
forall a. Num a => Integer -> a
fromInteger Integer
count
         in Double -> Either MathError Double
checked Double
mean Either MathError Double
-> Either MathError Double -> Either MathError Double
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Double -> Either MathError Double
checked Double
second Either MathError Double
-> Either MathError Streaming -> Either MathError Streaming
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Streaming -> Either MathError Streaming
forall a. a -> Either MathError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Integer -> Double -> Double -> Streaming
Streaming Integer
count Double
mean Double
second)

-- | Return count, mean and variance; True uses n-1 correction and requires at least two samples.
streamingResult :: Bool -> Streaming -> Either MathError (Integer, Double, Double)
streamingResult :: Bool -> Streaming -> Either MathError (Integer, Double, Double)
streamingResult Bool
correction (Streaming Integer
n Double
m Double
s) = do
    Bool -> String -> Either MathError ()
require (Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> if Bool
correction then Integer
1 else Integer
0) String
"insufficient streaming count"
    (Integer, Double, Double)
-> Either MathError (Integer, Double, Double)
forall a. a -> Either MathError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Integer
n, Double
m, Double
s Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Integer -> Double
forall a. Num a => Integer -> a
fromInteger (Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- if Bool
correction then Integer
1 else Integer
0))

{- | Histogram edges must increase strictly; intervals are [a,b), with the last
interval closed on its right. Points outside any axis are excluded. Flattening
is row-major, the last axis varies fastest. Weights need not sum to one.
-}
histogramND :: [[Double]] -> [[Double]] -> [Double] -> Either MathError [Double]
histogramND :: [[Double]] -> [[Double]] -> [Double] -> Either MathError [Double]
histogramND [[Double]]
edges [[Double]]
points [Double]
weights = do
    Bool -> String -> Either MathError ()
require
        ( Bool -> Bool
not ([[Double]] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [[Double]]
edges)
            Bool -> Bool -> Bool
&& ([Double] -> Bool) -> [[Double]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\[Double]
es -> [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
es Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
2 Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
es Bool -> Bool -> Bool
&& [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and ((Double -> Double -> Bool) -> [Double] -> [Double] -> [Bool]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
(<) [Double]
es (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
es))) [[Double]]
edges
        )
        String
"histogram edges invalid"
    Bool -> String -> Either MathError ()
require
        ( [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
points Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
weights
            Bool -> Bool -> Bool
&& ([Double] -> Bool) -> [[Double]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\[Double]
p -> [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
p Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
edges Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
p) [[Double]]
points
            Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Double
w -> Double -> Bool
finite Double
w Bool -> Bool -> Bool
&& Double
w Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0) [Double]
weights
        )
        String
"histogram observations/weights invalid"
    let dims :: [Int]
dims = ([Double] -> Int) -> [[Double]] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map ((Int -> Int -> Int
forall a. Num a => a -> a -> a
subtract Int
1) (Int -> Int) -> ([Double] -> Int) -> [Double] -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) [[Double]]
edges; size :: Integer
size = [Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ((Int -> Integer) -> [Int] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Integer
forall a. Integral a => a -> Integer
toInteger [Int]
dims)
    Bool -> String -> Either MathError ()
require (Integer
size Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
1000000) String
"histogram size budget exceeded"
    let index :: [Double] -> Maybe [Int]
index [Double]
p =
            [Maybe Int] -> Maybe [Int]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
forall (m :: * -> *) a. Monad m => [m a] -> m [a]
sequence
                [ if Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< [Double]
es [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0 Bool -> Bool -> Bool
|| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> [Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
es
                    then Maybe Int
forall a. Maybe a
Nothing
                    else Int -> Maybe Int
forall a. a -> Maybe a
Just (Int -> Int -> Int
forall a. Ord a => a -> a -> a
min ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
es Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2) ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((Double -> Bool) -> [Double] -> [Double]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
x) [Double]
es) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1))
                | ([Double]
es, Double
x) <- [[Double]] -> [Double] -> [([Double], Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [[Double]]
edges [Double]
p
                ]
        pairs :: [(Int, Double)]
pairs =
            [ ((Int -> (Int, Int) -> Int) -> Int -> [(Int, Int)] -> Int
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
Foldable.foldl' (\Int
a (Int
d, Int
i) -> Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
i) Int
0 ([Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int]
dims [Int]
ix), Double
w)
            | ([Double]
p, Double
w) <- [[Double]] -> [Double] -> [([Double], Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [[Double]]
points [Double]
weights
            , Just [Int]
ix <- [[Double] -> Maybe [Int]
index [Double]
p]
            ]
    (Double -> Either MathError Double)
-> [Double] -> Either MathError [Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM Double -> Either MathError Double
checked [[Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double
w | (Int
j, Double
w) <- [(Int, Double)]
pairs, Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i] | Int
i <- [Int
0 .. Integer -> Int
forall a. Num a => Integer -> a
fromInteger Integer
size Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]

validateProbabilities :: [Double] -> Either MathError ()
validateProbabilities :: [Double] -> Either MathError ()
validateProbabilities [Double]
ps =
    Bool -> String -> Either MathError ()
require
        (Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
ps) Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Double
p -> Double -> Bool
finite Double
p Bool -> Bool -> Bool
&& Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0) [Double]
ps Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
ps Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1e-10)
        String
"probabilities must be nonnegative and sum to one"

-- | Shannon entropy in nats for normalized nonnegative probabilities, with 0*log(0)=0.
entropy :: [Double] -> Either MathError Double
entropy :: [Double] -> Either MathError Double
entropy [Double]
ps = do [Double] -> Either MathError ()
validateProbabilities [Double]
ps; Double -> Either MathError Double
forall a. a -> Either MathError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (-[Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [if Double
p Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
0 else Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
p | Double
p <- [Double]
ps])

-- | Mutual information in nats from a rectangular normalized joint-probability table.
mutualInformation :: [[Double]] -> Either MathError Double
mutualInformation :: [[Double]] -> Either MathError Double
mutualInformation [[Double]]
joint = do
    Bool -> String -> Either MathError ()
require
        (Bool -> Bool
not ([[Double]] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [[Double]]
joint) Bool -> Bool -> Bool
&& Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null ([[Double]]
joint [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
0)) Bool -> Bool -> Bool
&& ([Double] -> Bool) -> [[Double]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([[Double]]
joint [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
0)) (Int -> Bool) -> ([Double] -> Int) -> [Double] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) [[Double]]
joint)
        String
"joint probability shape invalid"
    [Double] -> Either MathError ()
validateProbabilities ([[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Double]]
joint)
    let rows :: [Double]
rows = ([Double] -> Double) -> [[Double]] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [[Double]]
joint; cols :: [Double]
cols = ([Double] -> Double) -> [[Double]] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
transpose [[Double]]
joint)
    Double -> Either MathError Double
forall a. a -> Either MathError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure
        ( [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
            [ if Double
p Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
0 else Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log (Double
p Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ([Double]
rows [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
cols [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j))
            | (Int
i, [Double]
row) <- [Int] -> [[Double]] -> [(Int, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [[Double]]
joint
            , (Int
j, Double
p) <- [Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [Double]
row
            ]
        )

-- | KL(P||Q) in nats, including positive infinity when Q is zero on positive P mass.
klDivergence :: [Double] -> [Double] -> Either MathError Double
klDivergence :: [Double] -> [Double] -> Either MathError Double
klDivergence [Double]
p [Double]
q = do
    [Double] -> Either MathError ()
validateProbabilities [Double]
p
    [Double] -> Either MathError ()
validateProbabilities [Double]
q
    Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
p Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
q) String
"KL length mismatch"
    Double -> Either MathError Double
forall a. a -> Either MathError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [if Double
a Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
0 else if Double
b Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0 else Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log (Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
b) | (Double
a, Double
b) <- [Double] -> [Double] -> [(Double, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
p [Double]
q])

-- | Equal-mixture Jensen--Shannon divergence in nats, bounded by log 2.
jsDivergence :: [Double] -> [Double] -> Either MathError Double
jsDivergence :: [Double] -> [Double] -> Either MathError Double
jsDivergence [Double]
p [Double]
q = do
    let m :: [Double]
m = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
a Double
b -> (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) [Double]
p [Double]
q
    a <- [Double] -> [Double] -> Either MathError Double
klDivergence [Double]
p [Double]
m
    b <- klDivergence q m
    pure ((a + b) / 2)

-- | Standardize with population standard deviation; constant samples are rejected.
zScores :: [Double] -> Either MathError [Double]
zScores :: [Double] -> Either MathError [Double]
zScores [Double]
xs = do
    m <- [Double] -> Either MathError Double
arithmeticMean [Double]
xs
    s <- populationSD xs
    require (s > 0) "z-scores undefined for constant data"
    pure (map (\Double
x -> (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
m) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
s) xs)

-- | Unbiased sample standard deviation divided by sqrt n, requiring n>=2.
standardError :: [Double] -> Either MathError Double
standardError :: [Double] -> Either MathError Double
standardError [Double]
xs = (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs))) (Double -> Double)
-> Either MathError Double -> Either MathError Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [Double] -> Either MathError Double
sampleSD [Double]
xs

-- | Two-sided Student-t confidence interval for a normal population mean with unknown variance; confidence level in (0,1).
confidenceInterval :: Double -> [Double] -> Either MathError (Double, Double)
confidenceInterval :: Double -> [Double] -> Either MathError (Double, Double)
confidenceInterval Double
level [Double]
xs = do
    Bool -> String -> Either MathError ()
require (Double
level Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
level Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1) String
"confidence level outside (0,1)"
    m <- [Double] -> Either MathError Double
arithmeticMean [Double]
xs
    se <- standardError xs
    q <- P.quantile (P.StudentT (fromIntegral (length xs - 1))) ((1 + level) / 2)
    pure (m - q * se, m + q * se)

-- | Seeded nonparametric resampling with replacement. Returns caller statistic for each replicate and final RNG state.
bootstrap ::
    Int ->
    ([Double] -> Either MathError Double) ->
    [Double] ->
    R.RNG ->
    Either MathError ([Double], R.RNG)
bootstrap :: Int
-> ([Double] -> Either MathError Double)
-> [Double]
-> RNG
-> Either MathError ([Double], RNG)
bootstrap Int
count [Double] -> Either MathError Double
statistic [Double]
xs RNG
rng = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs
    Int
-> (RNG -> Either MathError (Double, RNG))
-> RNG
-> Either MathError ([Double], RNG)
forall a.
Int
-> (RNG -> Either MathError (a, RNG))
-> RNG
-> Either MathError ([a], RNG)
R.sampleN
        Int
count
        (\RNG
r -> do (ys, s) <- Int -> [Double] -> RNG -> Either MathError ([Double], RNG)
forall a. Int -> [a] -> RNG -> Either MathError ([a], RNG)
R.withReplacement ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs) [Double]
xs RNG
r; v <- statistic ys; pure (v, s))
        RNG
rng

-- | Leave-one-out bias-corrected statistic and jackknife standard error; requires n>=2.
jackknife :: ([Double] -> Either MathError Double) -> [Double] -> Either MathError (Double, Double)
jackknife :: ([Double] -> Either MathError Double)
-> [Double] -> Either MathError (Double, Double)
jackknife [Double] -> Either MathError Double
statistic [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
2 [Double]
xs
    theta <- [Double] -> Either MathError Double
statistic [Double]
xs
    ts <- mapM (\Int
i -> [Double] -> Either MathError Double
statistic (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
i [Double]
xs [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [Double]
xs)) [0 .. length xs - 1]
    m <- arithmeticMean ts
    let n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs)
     in pure (n * theta - (n - 1) * m, sqrt ((n - 1) / n * sum (map (\Double
t -> (Double
t 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
2 :: Int)) ts)))

-- | Test statistic, p value clamped to [0,1] after numerical checks, and applicable degrees of freedom.
data TestResult = TestResult
    { TestResult -> Double
testStatistic :: Double
    -- ^ Computed statistic for the selected hypothesis test.
    , TestResult -> Double
pValue :: Double
    -- ^ Reported tail probability, clamped to [0,1]; interpret using the test's documented alternative hypothesis.
    , TestResult -> [Double]
degreesOfFreedom :: [Double]
    -- ^ Distribution degrees-of-freedom parameters in test-specific order (for example numerator then denominator for F tests).
    }
    deriving (TestResult -> TestResult -> Bool
(TestResult -> TestResult -> Bool)
-> (TestResult -> TestResult -> Bool) -> Eq TestResult
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: TestResult -> TestResult -> Bool
== :: TestResult -> TestResult -> Bool
$c/= :: TestResult -> TestResult -> Bool
/= :: TestResult -> TestResult -> Bool
Eq, Int -> TestResult -> String -> String
[TestResult] -> String -> String
TestResult -> String
(Int -> TestResult -> String -> String)
-> (TestResult -> String)
-> ([TestResult] -> String -> String)
-> Show TestResult
forall a.
(Int -> a -> String -> String)
-> (a -> String) -> ([a] -> String -> String) -> Show a
$cshowsPrec :: Int -> TestResult -> String -> String
showsPrec :: Int -> TestResult -> String -> String
$cshow :: TestResult -> String
show :: TestResult -> String
$cshowList :: [TestResult] -> String -> String
showList :: [TestResult] -> String -> String
Show)

result :: Double -> Double -> [Double] -> Either MathError TestResult
result :: Double -> Double -> [Double] -> Either MathError TestResult
result Double
stat Double
p [Double]
df = do
    _ <- Double -> Either MathError Double
checked Double
stat
    require (finite p) "test p-value not finite"
    pure (TestResult stat (max 0 (min 1 p)) df)
tTail :: Double -> Double -> Either MathError Double
tTail :: Double -> Double -> Either MathError Double
tTail Double
df Double
t = (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
*) (Double -> Double)
-> Either MathError Double -> Either MathError Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Distribution -> Double -> Either MathError Double
P.survival (Double -> Distribution
P.StudentT Double
df) (Double -> Double
forall a. Num a => a -> a
abs Double
t)
chiTail :: Double -> Double -> Either MathError Double
chiTail :: Double -> Double -> Either MathError Double
chiTail Double
df Double
x = Distribution -> Double -> Either MathError Double
P.survival (Double -> Distribution
P.ChiSquare Double
df) Double
x

-- | Two-sided one-sample Student t test against a finite null mean; assumes independent normal observations.
oneSampleT :: Double -> [Double] -> Either MathError TestResult
oneSampleT :: Double -> [Double] -> Either MathError TestResult
oneSampleT Double
nullMean [Double]
xs = do
    m <- [Double] -> Either MathError Double
arithmeticMean [Double]
xs
    se <- standardError xs
    require (finite nullMean && se > 0) "t-test positive standard error required"
    let t = (Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
nullMean) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
se; df = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
     in do p <- tTail df t; result t p [df]

-- | Two-sided pooled t test when True, Welch unequal-variance test when False; each group requires n>=2 and positive combined variation.
twoSampleT :: Bool -> [Double] -> [Double] -> Either MathError TestResult
twoSampleT :: Bool -> [Double] -> [Double] -> Either MathError TestResult
twoSampleT Bool
equalVariance [Double]
xs [Double]
ys = do
    a <- [Double] -> Either MathError Double
arithmeticMean [Double]
xs
    b <- arithmeticMean ys
    v <- sampleVariance xs
    w <- sampleVariance ys
    let n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs)
        m = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys)
        s1 = Double
v Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
n
        s2 = Double
w Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
m
        se =
            if Bool
equalVariance
                then Double -> Double
forall a. Floating a => a -> a
sqrt (((Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
m))
                else Double -> Double
forall a. Floating a => a -> a
sqrt (Double
s1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
s2)
        df =
            if Bool
equalVariance
                then Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2
                else (Double
s1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
s2) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
s1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
s2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1))
        t = (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
se
    require (se > 0) "two-sample t requires positive variance"
    p <- tTail df t
    result t p [df]

-- | Two-sided one-sample t test on paired differences; matching group lengths required.
pairedT :: [Double] -> [Double] -> Either MathError TestResult
pairedT :: [Double] -> [Double] -> Either MathError TestResult
pairedT [Double]
xs [Double]
ys = do Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys) String
"paired lengths differ"; Double -> [Double] -> Either MathError TestResult
oneSampleT Double
0 ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Double]
xs [Double]
ys)

-- | Two-sided one-sample normal z test with known strictly positive population standard deviation.
zTest :: Double -> Double -> [Double] -> Either MathError TestResult
zTest :: Double -> Double -> [Double] -> Either MathError TestResult
zTest Double
nullMean Double
knownSD [Double]
xs = do
    m <- [Double] -> Either MathError Double
arithmeticMean [Double]
xs
    require (finite nullMean && finite knownSD && knownSD > 0) "z-test needs positive known SD"
    let z = (Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
nullMean) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
knownSD
     in do p <- (2 *) <$> P.survival (P.Normal 0 1) (abs z); result z p []

-- | Pearson goodness-of-fit chi-square statistic for observed and positive expected counts; first argument is number of fitted nuisance parameters.
chiSquareTest :: Int -> [Double] -> [Double] -> Either MathError TestResult
chiSquareTest :: Int -> [Double] -> [Double] -> Either MathError TestResult
chiSquareTest Int
fitted [Double]
xs [Double]
expected = do
    Int -> [Double] -> Either MathError ()
validate Int
2 [Double]
xs
    Int -> [Double] -> Either MathError ()
validate Int
2 [Double]
expected
    Bool -> String -> Either MathError ()
require
        ( [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
expected
            Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0) [Double]
xs
            Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) [Double]
expected
            Bool -> Bool -> Bool
&& Int
fitted Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0
            Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
fitted Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
        )
        String
"chi-square invalid counts/df"
    let stat :: Double
stat = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
a Double
e -> (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
e) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
e) [Double]
xs [Double]
expected)
        df :: Double
df = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
fitted Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
     in do p <- Double -> Double -> Either MathError Double
chiTail Double
df Double
stat; result stat p [df]

-- | Two-sided normal-population variance-ratio F test, requiring positive sample variances.
varianceFTest :: [Double] -> [Double] -> Either MathError TestResult
varianceFTest :: [Double] -> [Double] -> Either MathError TestResult
varianceFTest [Double]
xs [Double]
ys = do
    a <- [Double] -> Either MathError Double
sampleVariance [Double]
xs
    b <- sampleVariance ys
    require (a > 0 && b > 0) "variance F needs positive variances"
    let f = Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
b; df1 = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1); df2 = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
     in do p <- P.cdf (P.FisherF df1 df2) f; result f (2 * min p (1 - p)) [df1, df2]

validateGroups :: [[Double]] -> Either MathError ()
validateGroups :: [[Double]] -> Either MathError ()
validateGroups [[Double]]
gs = do Bool -> String -> Either MathError ()
require ([[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
gs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
2) String
"at least two groups required"; ([Double] -> Either MathError ())
-> [[Double]] -> Either MathError ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (Int -> [Double] -> Either MathError ()
validate Int
2) [[Double]]
gs

-- | One-way fixed-effects F test on independent normal groups with common variance; each group requires n>=2.
anova :: [[Double]] -> Either MathError TestResult
anova :: [[Double]] -> Either MathError TestResult
anova [[Double]]
gs = do
    [[Double]] -> Either MathError ()
validateGroups [[Double]]
gs
    overall <- [Double] -> Either MathError Double
arithmeticMean ([[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Double]]
gs)
    means <- mapM arithmeticMean gs
    let k = [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
gs
        n = [Int] -> Int
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (([Double] -> Int) -> [[Double]] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
gs)
        between = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
g) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
overall) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) | ([Double]
g, Double
m) <- [[Double]] -> [Double] -> [([Double], Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [[Double]]
gs [Double]
means]
        within = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [[Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Double
x -> (Double
x 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
2 :: Int)) [Double]
g) | ([Double]
g, Double
m) <- [[Double]] -> [Double] -> [([Double], Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [[Double]]
gs [Double]
means]
        df1 = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
        df2 = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
k)
        f = (Double
between Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
df1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
within Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
df2)
    require (within > 0) "ANOVA needs within-group variation"
    p <- P.survival (P.FisherF df1 df2) f
    result f p [df1, df2]

-- | Welch one-way heteroscedastic F test with Satterthwaite denominator degrees of freedom.
welchAnova :: [[Double]] -> Either MathError TestResult
welchAnova :: [[Double]] -> Either MathError TestResult
welchAnova [[Double]]
gs = do
    [[Double]] -> Either MathError ()
validateGroups [[Double]]
gs
    ms <- ([Double] -> Either MathError Double)
-> [[Double]] -> Either MathError [Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM [Double] -> Either MathError Double
arithmeticMean [[Double]]
gs
    vs <- mapM sampleVariance gs
    require (all (> 0) vs) "Welch ANOVA needs positive group variances"
    let ns = ([Double] -> Double) -> [[Double]] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Double) -> ([Double] -> Int) -> [Double] -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) [[Double]]
gs
        k = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
gs)
        ws = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Fractional a => a -> a -> a
(/) [Double]
ns [Double]
vs
        total = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
ws
        m = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) [Double]
ws [Double]
ms) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
total
        s = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [((Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
w Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
total) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) | (Double
w, Double
n) <- [Double] -> [Double] -> [(Double, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
ws [Double]
ns]
        f =
            ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
w Double
u -> Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
u 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
2 :: Int)) [Double]
ws [Double]
ms) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1))
                Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s)
        df = (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s)
    p <- P.survival (P.FisherF (k - 1) df) f
    result f p [k - 1, df]

tieCorrection :: [Double] -> Double
tieCorrection :: [Double] -> Double
tieCorrection [Double]
xs =
    let n :: Double
n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs)
     in Double
1
            Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
g) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
3 :: Int) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
g) | [Double]
g <- [Double] -> [[Double]]
forall a. Eq a => [a] -> [[a]]
group ([Double] -> [Double]
forall a. Ord a => [a] -> [a]
sort [Double]
xs)]
                Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
n)

-- | Tie-corrected Kruskal--Wallis H statistic with chi-square asymptotic calibration.
kruskalWallis :: [[Double]] -> Either MathError TestResult
kruskalWallis :: [[Double]] -> Either MathError TestResult
kruskalWallis [[Double]]
gs = do
    [[Double]] -> Either MathError ()
validateGroups [[Double]]
gs
    let xs :: [Double]
xs = [[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Double]]
gs
        rs :: [Double]
rs = [Double] -> [Double]
ranks [Double]
xs
        offsets :: [Int]
offsets = (Int -> Int -> Int) -> Int -> [Int] -> [Int]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl Int -> Int -> Int
forall a. Num a => a -> a -> a
(+) Int
0 (([Double] -> Int) -> [[Double]] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
gs)
        n :: Double
n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs)
        c :: Double
c = [Double] -> Double
tieCorrection [Double]
xs
        h :: Double
h =
            ( Double
12
                Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1))
                Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
                    [ ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
g) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
off [Double]
rs))) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
g)
                    | ([Double]
g, Int
off) <- [[Double]] -> [Int] -> [([Double], Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [[Double]]
gs [Int]
offsets
                    ]
                Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1)
            )
                Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
c
        df :: Double
df = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
gs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
    Bool -> String -> Either MathError ()
require (Double
c Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"Kruskal-Wallis requires nonconstant ranks"
    p <- Double -> Double -> Either MathError Double
chiTail Double
df Double
h
    result h p [df]

-- | Mann--Whitney U for the first sample; two-sided tie-corrected normal approximation with continuity correction.
mannWhitney :: [Double] -> [Double] -> Either MathError TestResult
mannWhitney :: [Double] -> [Double] -> Either MathError TestResult
mannWhitney [Double]
xs [Double]
ys = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
ys
    let n :: Double
n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs)
        m :: Double
m = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys)
        u :: Double
u = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs) ([Double] -> [Double]
ranks ([Double]
xs [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
ys))) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
        v :: Double
v = Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
12 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
tieCorrection ([Double]
xs [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
ys)
        z :: Double
z = (Double -> Double
forall a. Num a => a -> a
abs (Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
m Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt Double
v
     in do
            Bool -> String -> Either MathError ()
require (Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"Mann-Whitney constant ranks"
            p <- (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
*) (Double -> Double)
-> Either MathError Double -> Either MathError Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Distribution -> Double -> Either MathError Double
P.survival (Double -> Double -> Distribution
P.Normal Double
0 Double
1) (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 Double
z)
            result u p []

-- | Wilcoxon signed-rank test for differences, dropping zeros and using average absolute tied ranks and continuity-corrected normal calibration.
wilcoxon :: [Double] -> Either MathError TestResult
wilcoxon :: [Double] -> Either MathError TestResult
wilcoxon [Double]
differences = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
differences
    let xs :: [Double]
xs = (Double -> Bool) -> [Double] -> [Double]
forall a. (a -> Bool) -> [a] -> [a]
filter (Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
0) [Double]
differences
        n :: Double
n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs)
        rs :: [Double]
rs = [Double] -> [Double]
ranks ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs [Double]
xs)
        plus :: Double
plus = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double
r | (Double
x, Double
r) <- [Double] -> [Double] -> [(Double, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
xs [Double]
rs, Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0]
        mu :: Double
mu = Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4
        ties :: Double
ties =
            [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
                [Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
g) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
3 :: Int) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
g) | [Double]
g <- [Double] -> [[Double]]
forall a. Eq a => [a] -> [[a]]
group ([Double] -> [Double]
forall a. Ord a => [a] -> [a]
sort ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs [Double]
xs))]
        v :: Double
v = Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
24 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
ties Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
48
        z :: Double
z = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 (Double -> Double
forall a. Num a => a -> a
abs (Double
plus Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
mu) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt Double
v
    Bool -> String -> Either MathError ()
require (Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"Wilcoxon needs nonzero differences"
    p <- (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
*) (Double -> Double)
-> Either MathError Double -> Either MathError Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Distribution -> Double -> Either MathError Double
P.survival (Double -> Double -> Distribution
P.Normal Double
0 Double
1) Double
z
    result (min plus (n * (n + 1) / 2 - plus)) p []

-- | Two-sided KS statistic and asymptotic p value against a fully specified continuous distribution. A Jacobi theta transformation handles very small statistics. Fitting parameters from these observations invalidates this null calibration.
kolmogorovSmirnov :: P.Distribution -> [Double] -> Either MathError TestResult
kolmogorovSmirnov :: Distribution -> [Double] -> Either MathError TestResult
kolmogorovSmirnov Distribution
d [Double]
xs = do
    stat <- Distribution -> [Double] -> Either MathError Double
P.goodnessOfFit Distribution
d [Double]
xs
    let root = Double -> Double
forall a. Floating a => a -> a
sqrt (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs))
        z = (Double
root Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.12 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.11 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
root) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
stat
        -- Jacobi's theta transformation converges for small z, where a fixed
        -- alternating survival sum would need thousands or millions of terms.
        probability
            | Double
z Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 = Double
1
            | Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1.18 =
                Double
1
                    Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
sqrt (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi)
                        Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
z
                        Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double -> Double
forall a. Floating a => a -> a
exp (-Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int -> Int -> Int
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
8 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z)) | Int
k <- [Int
1 .. Int
20 :: Int]]
            | Bool
otherwise =
                Double
2
                    Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [(if Int -> Bool
forall a. Integral a => a -> Bool
odd Int
k then Double
1 else -Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp (-Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z) | Int
k <- [Int
1 .. Int
100 :: Int]]
    result stat probability []

{- | Royston's AS R94 coefficients and normal-tail approximation, n=3..5000,
uncensored samples. Numerical constants are published statistical data.
Reference: https://doi.org/10.2307/2986146
-}
shapiroWilk :: [Double] -> Either MathError TestResult
shapiroWilk :: [Double] -> Either MathError TestResult
shapiroWilk [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
3 [Double]
xs
    Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
5000) String
"Shapiro-Wilk supports n<=5000"
    let n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs
        an :: Double
an = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
        half :: Int
half = Int
n Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
        ys :: [Double]
ys = [Double] -> [Double]
forall a. Ord a => [a] -> [a]
sort [Double]
xs
        poly :: [a] -> a -> a
poly [a]
cs a
z = [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
z a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
i | (Int
i, a
c) <- [Int] -> [a] -> [(Int, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [a]
cs]
    v <- [Double] -> Either MathError Double
populationVariance [Double]
xs
    require (v > 0) "Shapiro-Wilk needs nonconstant observations"
    coeff <-
        if n == 3
            then pure [sqrt 0.5]
            else do
                zs <- mapM (\Int
i -> Double -> Either MathError Double
S.normalQuantile ((Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.375) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
an Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.25))) [1 .. half]
                let sumsq = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Double
z -> Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z) [Double]
zs)
                    scale = Double -> Double
forall a. Floating a => a -> a
sqrt Double
sumsq
                    r = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt Double
an
                    a1 = [Double] -> Double -> Double
forall {a}. Num a => [a] -> a -> a
poly [Double
0, Double
0.221157, -Double
0.147981, -Double
2.071190, Double
4.434685, -Double
2.706056] Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
zs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
scale
                    a2 = [Double] -> Double -> Double
forall {a}. Num a => [a] -> a -> a
poly [Double
0, Double
0.042981, -Double
0.293762, -Double
1.752461, Double
5.682633, -Double
3.582633] Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
zs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
scale
                    start = if Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
5 then Int
2 else Int
1
                    adjust =
                        if Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
5
                            then
                                Double -> Double
forall a. Floating a => a -> a
sqrt
                                    ((Double
sumsq Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
zs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
zs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
1) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a2))
                            else Double -> Double
forall a. Floating a => a -> a
sqrt ((Double
sumsq Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
zs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a1))
                pure (take start [a1, a2] ++ map (\Double
z -> -Double
z Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
adjust) (drop start zs))
    let numerator = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
ys [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
ys [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) | (Int
i, Double
a) <- [Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [Double]
coeff]
        w = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 (Double -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
1 (Double
numerator Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
numerator Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
an Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v)))
    p <-
        if n == 3
            then pure (1 - 6 / pi * acos (sqrt w))
            else do
                let y = Double -> Double
forall a. Floating a => a -> a
log (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1e-300 (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
w))
                    g = -Double
2.273 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.459 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
an
                    small = Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
11
                    transformed = if Bool
small then -Double -> Double
forall a. Floating a => a -> a
log (Double
g Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
y) else Double
y
                    mu =
                        if Bool
small
                            then [Double] -> Double -> Double
forall {a}. Num a => [a] -> a -> a
poly [Double
0.544, -Double
0.39978, Double
0.025054, -Double
0.0006714] Double
an
                            else [Double] -> Double -> Double
forall {a}. Num a => [a] -> a -> a
poly [-Double
1.5861, -Double
0.31082, -Double
0.083751, Double
0.0038915] (Double -> Double
forall a. Floating a => a -> a
log Double
an)
                    sd =
                        Double -> Double
forall a. Floating a => a -> a
exp
                            ( if Bool
small
                                then [Double] -> Double -> Double
forall {a}. Num a => [a] -> a -> a
poly [Double
1.3822, -Double
0.77857, Double
0.062767, -Double
0.0020322] Double
an
                                else [Double] -> Double -> Double
forall {a}. Num a => [a] -> a -> a
poly [-Double
0.4803, -Double
0.082676, Double
0.0030302] (Double -> Double
forall a. Floating a => a -> a
log Double
an)
                            )
                if small && y >= g then pure 1e-99 else P.survival (P.Normal mu sd) transformed
    result w p []

{- | Anderson--Darling normality test after fitting the sample mean and unbiased
sample standard deviation. The returned statistic is raw A²; Stephens' p-value
approximation uses A²*(1+0.75/n+2.25/n²). For adjusted A²>13, the published
approximation is outside its calibrated range and the extreme tail is
reported as zero instead of extrapolating its quadratic back upward.
-}
andersonDarling :: [Double] -> Either MathError TestResult
andersonDarling :: [Double] -> Either MathError TestResult
andersonDarling [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
5 [Double]
xs
    mean <- [Double] -> Either MathError Double
arithmeticMean [Double]
xs
    sd <- sampleSD xs
    let d = Double -> Double -> Distribution
P.Normal Double
mean Double
sd
    ps <- mapM (P.cdf d) (sort xs)
    let n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs
        an = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
        a =
            -Double
an
                Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
                    [ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall a. Floating a => a -> a
log (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1e-300 Double
probability) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
log (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1e-300 (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
ps [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1))))
                    | (Int
i, Double
probability) <- [Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [Double]
ps
                    ]
                    Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
an
        adjust = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.75 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
an Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2.25 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
an Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
an))
        p
            | Double
adjust Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0.2 = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
exp (-Double
13.436 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
101.14 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
adjust Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
223.73 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
adjust Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
adjust)
            | Double
adjust Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0.34 = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
exp (-Double
8.318 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
42.796 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
adjust Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
59.938 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
adjust Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
adjust)
            | Double
adjust Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0.6 = Double -> Double
forall a. Floating a => a -> a
exp (Double
0.9177 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
4.279 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
adjust Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1.38 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
adjust Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
adjust)
            | Double
adjust Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
13 = Double -> Double
forall a. Floating a => a -> a
exp (Double
1.2937 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
5.709 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
adjust Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.0186 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
adjust Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
adjust)
            | Bool
otherwise = Double
0
    result a p []

-- | One-way ANOVA on absolute deviations; True chooses group medians (Brown--Forsythe), False chooses means.
levene :: Bool -> [[Double]] -> Either MathError TestResult
levene :: Bool -> [[Double]] -> Either MathError TestResult
levene Bool
useMedian [[Double]]
gs = do
    [[Double]] -> Either MathError ()
validateGroups [[Double]]
gs
    centers <- ([Double] -> Either MathError Double)
-> [[Double]] -> Either MathError [Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM (if Bool
useMedian then [Double] -> Either MathError Double
median else [Double] -> Either MathError Double
arithmeticMean) [[Double]]
gs
    anova [map (\Double
x -> Double -> Double
forall a. Num a => a -> a
abs (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
c)) g | (g, c) <- zip gs centers]

-- | Bartlett homogeneity-of-variance test, appropriate for normal samples; positive group variances required.
bartlett :: [[Double]] -> Either MathError TestResult
bartlett :: [[Double]] -> Either MathError TestResult
bartlett [[Double]]
gs = do
    [[Double]] -> Either MathError ()
validateGroups [[Double]]
gs
    vs <- ([Double] -> Either MathError Double)
-> [[Double]] -> Either MathError [Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM [Double] -> Either MathError Double
sampleVariance [[Double]]
gs
    require (all (> 0) vs) "Bartlett requires positive group variances"
    let ds = ([Double] -> Double) -> [[Double]] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Double) -> ([Double] -> Int) -> [Double] -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Int -> Int
forall a. Num a => a -> a -> a
subtract Int
1 (Int -> Int) -> ([Double] -> Int) -> [Double] -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) [[Double]]
gs
        total = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
ds
        k = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
gs)
        pooled = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) [Double]
ds [Double]
vs) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
total
        stat =
            (Double
total Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
pooled Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
d Double
v -> Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
v) [Double]
ds [Double]
vs))
                Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/) [Double]
ds) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
total) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)))
     in do p <- chiTail (k - 1) stat; result stat p [k - 1]

-- | Two-sided Fisher exact test on nonnegative 2-by-2 counts (a,b;c,d). Enumerates tables with mass no larger than the observed table.
fisherExact :: Int -> Int -> Int -> Int -> Either MathError TestResult
fisherExact :: Int -> Int -> Int -> Int -> Either MathError TestResult
fisherExact Int
a Int
b Int
c Int
d = do
    Bool -> String -> Either MathError ()
require
        ((Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0) [Int
a, Int
b, Int
c, Int
d] Bool -> Bool -> Bool
&& Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
b Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
b Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
100000)
        String
"Fisher counts invalid"
    let dist :: Distribution
dist = Int -> Int -> Int -> Distribution
P.Hypergeometric (Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
b Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
d) (Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c) (Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
b)
        lo :: Int
lo = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
0 (Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
d)
        hi :: Int
hi = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min (Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
b) (Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c)
     in do
            p0 <- Distribution -> Double -> Either MathError Double
P.density Distribution
dist (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
a)
            ps <- mapM (P.density dist . fromIntegral) [lo .. hi]
            let odds =
                    if Int
b Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
c Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
                        then if Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then Double
0 else Double
1e300
                        else Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
d Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
c)
             in result odds (sum (filter (<= p0 * (1 + 1e-12)) ps)) []

-- | Seeded two-sided random-label permutation p value for a caller statistic, with the (b+1)/(B+1) correction.
permutationTest ::
    Int ->
    ([Double] -> [Double] -> Either MathError Double) ->
    [Double] ->
    [Double] ->
    R.RNG ->
    Either MathError (TestResult, R.RNG)
permutationTest :: Int
-> ([Double] -> [Double] -> Either MathError Double)
-> [Double]
-> [Double]
-> RNG
-> Either MathError (TestResult, RNG)
permutationTest Int
count [Double] -> [Double] -> Either MathError Double
statistic [Double]
xs [Double]
ys RNG
rng = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
ys
    Bool -> String -> Either MathError ()
require (Int
count Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
count Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
1000000) String
"permutation count outside budget"
    observed <- [Double] -> [Double] -> Either MathError Double
statistic [Double]
xs [Double]
ys
    (values, s) <-
        R.sampleN
            count
            ( \RNG
r -> do
                (z, t) <- [Double] -> RNG -> Either MathError ([Double], RNG)
forall a. [a] -> RNG -> Either MathError ([a], RNG)
R.shuffle ([Double]
xs [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
ys) RNG
r
                v <- statistic (take (length xs) z) (drop (length xs) z)
                pure (v, t)
            )
            rng
    res <-
        result
            observed
            (fromIntegral (1 + length (filter (\Double
v -> Double -> Double
forall a. Num a => a -> a
abs Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double -> Double
forall a. Num a => a -> a
abs Double
observed) values)) / fromIntegral (count + 1))
            []
    pure (res, s)

-- | Bonferroni-adjusted p values in original order, capped at one.
bonferroni :: [Double] -> Either MathError [Double]
bonferroni :: [Double] -> Either MathError [Double]
bonferroni [Double]
ps = do
    Bool -> String -> Either MathError ()
require ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Double
p -> Double -> Bool
finite Double
p Bool -> Bool -> Bool
&& Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 Bool -> Bool -> Bool
&& Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1) [Double]
ps) String
"invalid p-values"
    [Double] -> Either MathError [Double]
forall a. a -> Either MathError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
1 (Double -> Double) -> (Double -> Double) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ps))) [Double]
ps)

-- | Monotone Benjamini--Hochberg FDR-adjusted p values in original order.
benjaminiHochberg :: [Double] -> Either MathError [Double]
benjaminiHochberg :: [Double] -> Either MathError [Double]
benjaminiHochberg [Double]
ps = do
    _ <- [Double] -> Either MathError [Double]
bonferroni [Double]
ps
    let ordered = ((Int, Double) -> Double) -> [(Int, Double)] -> [(Int, Double)]
forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn (Int, Double) -> Double
forall a b. (a, b) -> b
snd ([Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [Double]
ps)
        n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ps
        adjusted =
            [(Int, Double)] -> [(Int, Double)]
forall a. [a] -> [a]
reverse
                ( (Double, [(Int, Double)]) -> [(Int, Double)]
forall a b. (a, b) -> b
snd
                    ( ((Double, [(Int, Double)])
 -> (Int, (Int, Double)) -> (Double, [(Int, Double)]))
-> (Double, [(Int, Double)])
-> [(Int, (Int, Double))]
-> (Double, [(Int, Double)])
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
Foldable.foldl'
                        ( \(Double
previous, [(Int, Double)]
acc) (Int
rank, (Int
i, Double
p)) -> let q :: Double
q = Double -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
previous (Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
rank) in (Double
q, (Int
i, Double
q) (Int, Double) -> [(Int, Double)] -> [(Int, Double)]
forall a. a -> [a] -> [a]
: [(Int, Double)]
acc)
                        )
                        (Double
1, [])
                        ([(Int, (Int, Double))] -> [(Int, (Int, Double))]
forall a. [a] -> [a]
reverse ([Int] -> [(Int, Double)] -> [(Int, (Int, Double))]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
1 :: Int ..] [(Int, Double)]
ordered))
                    )
                )
    pure (map snd (sortOn fst adjusted))

-- | Difference of sample means divided by pooled unbiased standard deviation.
cohensD :: [Double] -> [Double] -> Either MathError Double
cohensD :: [Double] -> [Double] -> Either MathError Double
cohensD [Double]
xs [Double]
ys = do
    a <- [Double] -> Either MathError Double
arithmeticMean [Double]
xs
    b <- arithmeticMean ys
    v <- sampleVariance xs
    w <- sampleVariance ys
    let pooled =
            ( (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w)
                Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
+ [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2)
            )
     in require (pooled > 0) "Cohen d needs nonzero pooled variance" >> checked ((a - b) / sqrt pooled)

-- | One-way ANOVA explained sum of squares divided by total sum of squares.
etaSquared :: [[Double]] -> Either MathError Double
etaSquared :: [[Double]] -> Either MathError Double
etaSquared [[Double]]
gs = do
    [[Double]] -> Either MathError ()
validateGroups [[Double]]
gs
    m <- [Double] -> Either MathError Double
arithmeticMean ([[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Double]]
gs)
    ms <- mapM arithmeticMean gs
    let between = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
g) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
u 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
2 :: Int) | ([Double]
g, Double
u) <- [[Double]] -> [Double] -> [([Double], Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [[Double]]
gs [Double]
ms]
        total = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Double
x -> (Double
x 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
2 :: Int)) ([[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Double]]
gs))
     in require (total > 0) "eta squared needs nonconstant data" >> checked (between / total)

-- | Exact power under a known-variance two-sided one-sample normal z-test.
normalTestPower :: Double -> Double -> Int -> Either MathError Double
normalTestPower :: Double -> Double -> Int -> Either MathError Double
normalTestPower Double
alpha Double
standardizedEffect Int
n = do
    Bool -> String -> Either MathError ()
require
        (Double
alpha Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
alpha Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
standardizedEffect Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0)
        String
"power invalid alpha/effect/n"
    q <- Double -> Either MathError Double
S.normalQuantile (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
alpha Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)
    let shift = Double
standardizedEffect Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)
     in do a <- S.normalCDF (-q - shift); b <- P.survival (P.Normal 0 1) (q - shift); pure (a + b)

-- | Smallest integer n whose exact known-variance two-sided normal-test power meets the target; hard budget n<=100000000.
normalSampleSize :: Double -> Double -> Double -> Either MathError Int
normalSampleSize :: Double -> Double -> Double -> Either MathError Int
normalSampleSize Double
alpha Double
power Double
effect = do
    Bool -> String -> Either MathError ()
require
        (Double
alpha Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
alpha Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1 Bool -> Bool -> Bool
&& Double
power Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
alpha Bool -> Bool -> Bool
&& Double
power Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
effect Bool -> Bool -> Bool
&& Double
effect Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
0)
        String
"sample size invalid alpha/power/effect"
    let go :: Int -> Int -> Either MathError Int
go Int
lo Int
hi
            | Int
hi Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
100000000 = do
                p <- Double -> Double -> Int -> Either MathError Double
normalTestPower Double
alpha Double
effect Int
100000000
                if p >= power then search lo 100000000 else Left (NonConvergence "sample size budget" 100000000)
            | Bool
otherwise = do p <- Double -> Double -> Int -> Either MathError Double
normalTestPower Double
alpha Double
effect Int
hi; if p >= power then search lo hi else go hi (hi * 2)
        search :: Int -> Int -> Either MathError Int
search Int
lo Int
hi
            | Int
hi Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
lo Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
1 = Int -> Either MathError Int
forall a. a -> Either MathError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure Int
hi
            | Bool
otherwise =
                let m :: Int
m = Int
lo Int -> Int -> Int
forall a. Num a => a -> a -> a
+ (Int
hi Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
lo) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
                 in do p <- Double -> Double -> Int -> Either MathError Double
normalTestPower Double
alpha Double
effect Int
m; if p >= power then search lo m else search m hi
    p1 <- Double -> Double -> Int -> Either MathError Double
normalTestPower Double
alpha Double
effect Int
1
    if p1 >= power then pure 1 else go 1 2

-- | Zero-based indices outside [Q1-k*IQR,Q3+k*IQR], with positive multiplier k and type-7 quartiles.
tukeyOutliers :: Double -> [Double] -> Either MathError [Int]
tukeyOutliers :: Double -> [Double] -> Either MathError [Int]
tukeyOutliers Double
multiplier [Double]
xs = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
multiplier Bool -> Bool -> Bool
&& Double
multiplier Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"Tukey multiplier must be positive"
    q1 <- Double -> [Double] -> Either MathError Double
quantile Double
0.25 [Double]
xs
    q3 <- quantile 0.75 xs
    let iqr = Double
q3 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
q1
     in pure [i | (i, x) <- zip [0 ..] xs, x < q1 - multiplier * iqr || x > q3 + multiplier * iqr]

-- | Two-sided single-outlier Grubbs test, returning candidate index, G statistic and rejection decision at alpha; normality assumption required.
grubbs :: Double -> [Double] -> Either MathError (Int, Double, Bool)
grubbs :: Double -> [Double] -> Either MathError (Int, Double, Bool)
grubbs Double
alpha [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
3 [Double]
xs
    Bool -> String -> Either MathError ()
require (Double
alpha Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
alpha Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1) String
"Grubbs alpha outside (0,1)"
    m <- [Double] -> Either MathError Double
arithmeticMean [Double]
xs
    s <- sampleSD xs
    require (s > 0) "Grubbs needs nonconstant data"
    let n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs)
        distances = (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Double
x -> Double -> Double
forall a. Num a => a -> a
abs (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
m) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
s) [Double]
xs
        g = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
distances
        i = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((Double -> Bool) -> [Double] -> [Double]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
g) [Double]
distances)
     in do
            t <- P.quantile (P.StudentT (n - 2)) (1 - alpha / (2 * n))
            let critical = (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t)) in pure (i, g, g > critical)

-- | Keep True-mask observations, then validate retained values. Excluded NaN/infinity do not contaminate results.
masked :: [Bool] -> [Double] -> Either MathError [Double]
masked :: [Bool] -> [Double] -> Either MathError [Double]
masked [Bool]
mask [Double]
xs = do
    Bool -> String -> Either MathError ()
require ([Bool] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Bool]
mask Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs) String
"mask length mismatch"
    let ys :: [Double]
ys = [Double
x | (Bool
keep, Double
x) <- [Bool] -> [Double] -> [(Bool, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Bool]
mask [Double]
xs, Bool
keep] in Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
ys Either MathError ()
-> Either MathError [Double] -> Either MathError [Double]
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> [Double] -> Either MathError [Double]
forall a. a -> Either MathError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure [Double]
ys

-- | Huber M location by iteratively reweighted means, scaled by 1.4826*MAD; zero MAD returns the median. Positive tuning and finite data required.
huberLocation :: Double -> [Double] -> Either MathError Double
huberLocation :: Double -> [Double] -> Either MathError Double
huberLocation Double
tuning [Double]
xs = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
tuning Bool -> Bool -> Bool
&& Double
tuning Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"Huber tuning positive required"
    m <- [Double] -> Either MathError Double
median [Double]
xs
    mad <- medianAbsoluteDeviation xs
    let scale = Double
1.4826 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
mad
        go :: Int -> Double -> Either MathError Double
        go Int
0 Double
_ = MathError -> Either MathError Double
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"Huber location" Int
1000)
        go Int
n Double
u = do
            let ws :: [Double]
ws = (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Double
x -> let r :: Double
r = Double -> Double
forall a. Num a => a -> a
abs ((Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
u) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
scale) in if Double
r Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
tuning then Double
1 else Double
tuning Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
r) [Double]
xs
            v <- [Double] -> [Double] -> Either MathError Double
weightedMean [Double]
xs [Double]
ws
            if abs (v - u) <= 1e-10 * max 1 (abs u) then pure v else go (n - 1) v
    if mad == 0 then pure m else go 1000 m