module Amby.Numeric
(
contDistrDomain
, contDistrRange
, linspace
, arange
, random
, scoreAtPercentile
, interquartileRange
, freedmanDiaconisBins
) where
import Data.Vector.Generic ((!))
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector as V
import qualified Data.Vector.Algorithms.Intro as V
import Statistics.Distribution
import System.Random.MWC (withSystemRandom, asGenST)
contDistrDomain :: (ContDistr d) => d -> Int -> U.Vector Double
contDistrDomain d num = linspace (quantile d 0.0001) (quantile d 0.9999) num
contDistrRange :: (ContDistr d) => d -> U.Vector Double -> U.Vector Double
contDistrRange d x = U.map (density d) x
linspace :: Double -> Double -> Int -> U.Vector Double
linspace start stop num
| num < 0 = error ("Number of samples, " ++ show num ++ ", must be non-negative.")
| num == 0 || num == 1 = addStart $ U.generate num ((* delta) . fromIntegral)
| otherwise = addStart $ U.generate num ((* step) . fromIntegral)
where
delta = stop start
step = delta / fromIntegral (num 1)
addStart = U.map (realToFrac . (+ start))
arange :: Double -> Double -> Double -> U.Vector Double
arange start stop step = U.fromList [start,(start + step)..stop]
random :: (ContGen d) => d -> Int -> IO (U.Vector Double)
random d n = withSystemRandom . asGenST $ \gen ->
U.replicateM n $ genContVar d gen
scoreAtPercentile :: (G.Vector v Double) => v Double -> Double -> Double
scoreAtPercentile xs p
| n == 0 = modErr "scoreAtPercentile" "Percentile sample size is 0"
| p < 0 || p > 100 = modErr "scoreAtPercentile" "Percentile must be in [0, 100]"
| otherwise = (getScore . G.modify V.sort) xs
where
n = G.length xs
idx = (p / 100) * fromIntegral (n 1)
i = floor idx
j = ceiling idx
isInt :: Double -> Bool
isInt a = fromIntegral (round a :: Int) == a
interop :: (G.Vector v Double) => v Double -> Double
interop vec =
(vec ! i) * (fromIntegral j idx) +
(vec ! j) * (idx fromIntegral i)
getScore vec = if isInt idx
then vec ! i
else interop vec
interquartileRange :: (G.Vector v Double) => v Double -> Double
interquartileRange xs = scoreAtPercentile xs 75 scoreAtPercentile xs 25
freedmanDiaconisBins :: (G.Vector v Double) => v Double -> Int
freedmanDiaconisBins xs =
if h == 0
then floor $ sqrt $ (fromIntegral n :: Double)
else ceiling $ (G.maximum xs G.minimum xs) / h
where
n = G.length xs
h = 2 * interquartileRange xs / fromIntegral n ** (1 / 3)
modErr :: String -> String -> a
modErr f err = error
$ showString "Amby.Numeric."
$ showString f err