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
case f z `compare` 0 of
GT -> normal (zdz) 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
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)
kDecomp n d = (k, h)
where
dn = d * fromIntegral n
k = ceiling dn
h = (fromIntegral k dn)
kScale n = foldl f (1,0) [1..n]
where
f (s,e) i = shift (s * fromIntegral i / fromIntegral n, e)
kCdfMat n d = (matrix m m hMatCell, k)
where
spine n = (1 h^n) / factorial n
hMatCell r 0 | r == m1 = (1 2 * h ^ m + max 0 ((2 * h 1)^m)) / factorial m
hMatCell r 0 = spine (fromIntegral r + 1)
hMatCell r c | r == m1 = spine (m fromIntegral c)
hMatCell r c
| i j + 1 >= 0 = 1 / factorial (fromIntegral (ij+1))
| otherwise = 0
where i = r + 1; j = c + 1
(k, h) = kDecomp n d
m = 2 * k 1
kCdf n d
| d <= 0 = 0
| otherwise = multAndUnshift (kScale n) (indexM hN (k1) (k1), expShift)
where
(hN, expShift) = mPower hMat n
(hMat, k) = kCdfMat n d
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
mPower m 1 = (m, 0)
mPower m n
| even n = square (mPower m (n `div` 2))
| otherwise = m `mult` mPower m (n1)
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)