module Amby.Numeric
(
contDistrDomain
, contDistrRange
, linspace
, arange
, random
, scoreAtPercentile
, interquartileRange
, freedmanDiaconisBins
) where
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector as V
import Statistics.Distribution
import qualified Statistics.Quantile as Stats
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 -> Int -> Double
scoreAtPercentile xs p
| G.length xs == 0 =
modErr "scoreAtPercentile" "Cannot find percentile of empyt list"
| otherwise = Stats.weightedAvg p 100 xs
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
$ showString ": " err