ks-test (empty) → 0.1
raw patch · 5 files changed
+303/−0 lines, 5 filesdep +basedep +gammadep +random-fusetup-changed
Dependencies added: base, gamma, random-fu, roots, vector
Files
- Setup.lhs +5/−0
- ks-test.cabal +32/−0
- src/Data/Random/Distribution/Kolmogorov.hs +138/−0
- src/Math/Statistics/KSTest.hs +92/−0
- src/Numeric/LinearAlgebra.hs +36/−0
+ Setup.lhs view
@@ -0,0 +1,5 @@+#!/usr/bin/env runhaskell++> import Distribution.Simple+> main = defaultMain+
+ ks-test.cabal view
@@ -0,0 +1,32 @@+name: ks-test+version: 0.1+stability: provisional++cabal-version: >= 1.6+build-type: Simple++author: James Cook <mokus@deepbondi.net>+maintainer: James Cook <mokus@deepbondi.net>+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
+ src/Data/Random/Distribution/Kolmogorov.hs view
@@ -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+
+ src/Math/Statistics/KSTest.hs view
@@ -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+
+ src/Numeric/LinearAlgebra.hs view
@@ -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]]++