diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,5 @@ +#!/usr/bin/env runhaskell + +> import Distribution.Simple +> main = defaultMain + diff --git a/ks-test.cabal b/ks-test.cabal new file mode 100644 --- /dev/null +++ b/ks-test.cabal @@ -0,0 +1,32 @@ +name: ks-test +version: 0.1 +stability: provisional + +cabal-version: >= 1.6 +build-type: Simple + +author: James Cook +maintainer: James Cook +license: PublicDomain + +category: Math, Numerical +synopsis: Kolmogorov distribution and Kolmogorov-Smirnov test. +description: Kolmogorov distribution and Kolmogorov-Smirnov test. + +tested-with: GHC == 6.10.4, + GHC == 6.12.1, GHC == 6.12.3 + +source-repository head + type: darcs + location: http://code.haskell.org/~mokus/ks-test + +Library + hs-source-dirs: src + exposed-modules: Data.Random.Distribution.Kolmogorov + Math.Statistics.KSTest + other-modules: Numeric.LinearAlgebra + build-depends: base >= 3 && <5, + gamma, + random-fu >= 0.1 && < 0.3, + roots, + vector diff --git a/src/Data/Random/Distribution/Kolmogorov.hs b/src/Data/Random/Distribution/Kolmogorov.hs new file mode 100644 --- /dev/null +++ b/src/Data/Random/Distribution/Kolmogorov.hs @@ -0,0 +1,138 @@ +{-# LANGUAGE + MultiParamTypeClasses, + FlexibleInstances, + FlexibleContexts + #-} +-- |CDF of Kolmogorov's D-statistic, parameterized by sample size +module Data.Random.Distribution.Kolmogorov (D(..), kCdf, kCdfQuick)where + +import Numeric.LinearAlgebra +import Math.Gamma +import Math.Root.Bracket +import Math.Root.Finder +import Math.Root.Finder.Brent + +import Data.Random + +data D a = D Int + deriving (Eq, Show) + +instance Distribution D Double where + rvar (D n) = do + u <- stdUniform + let f x = kCdfQuick n x - u + eps = encodeFloat 1 (1 - floatDigits eps) + (x0,x1) = last (bracket f 0.01 0.1) + + case findRoot f x0 x1 eps :: Either (Brent Double Double) (Brent Double Double) of + Left stopped -> fail ("Failed to find root, final state was: " ++ show stopped) + Right root -> do + let z = estimateRoot root + dz = 0.5 * estimateError root + + -- Very slight "blur" to account for the error in the root. + -- Center the blur above or below the reported root depending + -- on whether the root is above or below the zero. + case f z `compare` 0 of + GT -> normal (z-dz) dz + EQ -> normal z dz + LT -> normal (z+dz) dz + + + +instance CDF D Double where + cdf (D n) = kCdfQuick n + +shiftPointU :: Num a => a +shiftPointU = 2 ^ 450 +shiftPointL :: Fractional a => a +shiftPointL = 2 ^^ (-450) + +shiftDist :: Integral a => a +shiftDist = 450 + +shiftBase :: Num a => a +shiftBase = 2 + +-- Shifting and unshifting operations: These manage the external exponent +-- for a value. +shift x = shiftBy id (*) x +shiftBy f (*) (x,e) + | fx < shiftPointL = (shiftPointU * x, e - shiftDist) + | fx > shiftPointU = (shiftPointL * x, e + shiftDist) + | otherwise = (x,e) + where fx = f x + +unshift (x,0) = x +unshift (x,e) + | e >= shiftDist = unshift (shiftPointU * x, e - shiftDist) + | e <= negate shiftDist = unshift (shiftPointL * x, e + shiftDist) + | otherwise = x * shiftBase ** realToFrac e + +multAndUnshift (a,aE) (b,bE) = unshift (a*b, aE+bE) + +-- Decompose d into k and h as described in the paper +kDecomp n d = (k, h) + where + dn = d * fromIntegral n + k = ceiling dn + h = (fromIntegral k - dn) + +-- Compute the scale factor n!/n^n in "shifted" form +kScale n = foldl f (1,0) [1..n] + where + f (s,e) i = shift (s * fromIntegral i / fromIntegral n, e) + +-- H matrix used in kCdf and k value +kCdfMat n d = (matrix m m hMatCell, k) + where + spine n = (1 - h^n) / factorial n + + hMatCell r 0 | r == m-1 = (1 - 2 * h ^ m + max 0 ((2 * h - 1)^m)) / factorial m + hMatCell r 0 = spine (fromIntegral r + 1) + hMatCell r c | r == m-1 = spine (m - fromIntegral c) + hMatCell r c + | i - j + 1 >= 0 = 1 / factorial (fromIntegral (i-j+1)) + | otherwise = 0 + where i = r + 1; j = c + 1 + + (k, h) = kDecomp n d + m = 2 * k - 1 + +-- CDF of Kolmogorov's D-statistic for a given sample size n +kCdf n d + | d <= 0 = 0 + | otherwise = multAndUnshift (kScale n) (indexM hN (k-1) (k-1), expShift) + where + (hN, expShift) = mPower hMat n + (hMat, k) = kCdfMat n d + +-- |'kCdf' with Marsaglia's long-computation shortcut approximation. +-- Accurate to about 7 decimal places in the right tail of the distribution. +kCdfQuick n d + | s > 7.24 || (s > 3.76 && n > 99) + = 1 - 2 * exp (negate (2.000071 + 0.331/sqrt n' + 1.409/n') * s) + | otherwise = kCdf n d + where s = d*d*n'; n' = fromIntegral n + +-- |Matrix power to Int exponent with explicitly managed exponent. Returns +-- a multiple of m^n, along with the natural logarithm of the factor. +-- +-- That is, if @(mn, logS) = mPower m n@ then m^n = (e^logS) * mn +mPower m 1 = (m, 0) +mPower m n + | even n = square (mPower m (n `div` 2)) + | otherwise = m `mult` mPower m (n-1) + where + k = min (matRows m) (matCols m) `div` 2 + + square (m, e) = m `mult` (m,e + e) + mult m1 (m2, e2) = shiftBy (\m -> indexM m k k) scale (multiply m1 m2, e2) + + +-- Analytic limiting form (as n -> ∞) +-- kCdfLim d = sqrt (2 * pi) * recip_x * sum [exp (negate ((2 * fromIntegral i - 1)^2) * pi_2 * 0.125 * recip_x*recip_x) | i <- [1..]] +-- where +-- recip_x = recip x +-- pi_2 = pi*pi + diff --git a/src/Math/Statistics/KSTest.hs b/src/Math/Statistics/KSTest.hs new file mode 100644 --- /dev/null +++ b/src/Math/Statistics/KSTest.hs @@ -0,0 +1,92 @@ +{-# LANGUAGE + MultiParamTypeClasses, + FlexibleContexts, FlexibleInstances, + GADTs, + ParallelListComp + #-} +module Math.Statistics.KSTest (ks, ksTest, KS(..)) where + +import Control.Monad +import Data.Random +import Data.Random.Distribution.Kolmogorov (kCdf, kCdfQuick) +import Data.List (sort) + +-- |Kolmogorov-Smirnov statistic for a set of data relative to a (continuous) +-- distribution with the given CDF. Returns 3 common forms of the statistic: +-- (K+, K-, D), where K+ and K- are Smirnov's one-sided forms as presented in +-- Knuth's Semi-Numerical Algorithms (TAOCP, vol. 2) and D is Kolmogorov's +-- undirected version. +-- +-- In particular, +-- +-- * K+ = sup(\x -> F_n(x) - F(x)) +-- * K- = sup(\x -> F(x) - F_n(x)) +-- * D = sup(\x -> abs(F_n(x) - F(x))) +-- +ks f n xs = (kPlus, kMinus, d) + where + sqrt_n = sqrt (fromIntegral n) + sorted = sort (map f xs) + + kPlus = maximum [ j//n - fx | j <- [1..] | fx <- sorted] + kMinus = maximum [ fx - j//n | j <- [0..] | fx <- sorted] + d = max kPlus kMinus + + infixl 7 // + x//y = fromIntegral x / fromIntegral y + +-- | @ksTest cdf xs@ +-- Computes the probability of a random data set (of the same size as xs) +-- drawn from a continuous distribution with the given CDF having the same +-- Kolmogorov statistic as xs. +-- +-- The statistic is the greatest absolute deviation of the empirical CDF of +-- XS from the assumed CDF @cdf@. +-- +-- If the data were, in fact, drawn from a distribution with the given CDF, +-- then the resulting p-value should be uniformly distributed over (0,1]. +ksTest f xs = 1 - kCdf n d + where + n = length xs + + (kPlus, kMinus, d) = ks f n xs + +-- |'KS' distribution: not really a standard mathematical concept, but still +-- a nice conceptual shift. @KS n d@ is the distribution of a random +-- variable constructed as a list of @n@ independent random variables of +-- distribution @d@. +-- +-- The corresponding 'CDF' instance implements the K-S test for such lists. +-- For example, if @xs@ is a list of length 100 believed to contain Beta(2,5) +-- variates, then @cdf (KS 100 (Beta 2 5))@ is the K-S test for that distribution. +-- (Note that if @length xs@ is not 100, then the result will be 0 because +-- such lists cannot arise from that 'KS' distribution. Somewhat arbitrarily, +-- all lists of \"impossible\" length are grouped at the bottom of the ordering +-- encoded by the 'CDF' instance.) +-- +-- The 'KS' test can easily be applied recursively. +-- For example, if @d@ is a 'Distribution' of interest and you have 100 trials +-- each with 100 data points, you can test it by calling @cdf (KS 100 (KS 100 d))@. +data KS d a where + KS :: !Int -> !(d a) -> KS d [a] + +instance Eq (d a) => Eq (KS d [a]) where + KS n1 d1 == KS n2 d2 = n1 == n2 && d1 == d2 +instance Show (d a) => Show (KS d [a]) where + showsPrec p (KS n d) = showParen (p > 10) + ( showString "KS " + . showsPrec 11 n + . showChar ' ' + . showsPrec 11 d + ) + +instance Distribution d a => Distribution (KS d) [a] where + rvar (KS n d) = replicateM n (rvar d) + +instance CDF d a => CDF (KS d) [a] where + cdf (KS n dist) xs + | length (take (n+1) xs) == n = 1 - kCdf n d + | otherwise = 0 + where + (kPlus, kMinus, d) = ks (cdf dist) n xs + diff --git a/src/Numeric/LinearAlgebra.hs b/src/Numeric/LinearAlgebra.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/LinearAlgebra.hs @@ -0,0 +1,36 @@ +-- |Minimal linear algebra support for implementing KS statistic CDF +module Numeric.LinearAlgebra where + +import Data.Vector.Unboxed (Vector, Unbox, generate, (!), map) + +data Matrix a = Matrix + { matRows :: !Int + , matCols :: !Int + , content :: !(Vector a) + } + +matrix :: Unbox a => Int -> Int -> (Int -> Int -> a) -> Matrix a +matrix r c f = Matrix r c (generate n (uncurry f . idx)) + where + n = r*c + idx i = i `divMod` r + +indexM :: Unbox a => Matrix a -> Int -> Int -> a +indexM m@(Matrix r c v) i j + | i < 0 = error "indexM: i < 0" + | j < 0 = error "indexM: j < 0" + | i >= r = error "indexM: i >= r" + | j >= c = error "indexM: j >= c" + | otherwise = unsafeIndexM m i j + +unsafeIndexM (Matrix r c v) i j = v ! (i * r + j) + +scale :: (Unbox a, Num a) => a -> Matrix a -> Matrix a +scale k (Matrix r c v) = Matrix r c (Data.Vector.Unboxed.map (*k) v) + +multiply :: (Unbox a, Num a) => Matrix a -> Matrix a -> Matrix a +multiply m1@(Matrix r1 c1 _) m2@(Matrix r2 c2 _) + | c1 /= r2 = error "multiply: incompatible matrix sizes" + | otherwise = matrix r1 c2 $ \i j -> sum [unsafeIndexM m1 i k * unsafeIndexM m2 k j | k <- [0..c1-1]] + +