statistics 0.6.0.0 → 0.6.0.1
raw patch · 3 files changed
+203/−48 lines, 3 filesdep ~mwc-randomdep ~primitivedep ~vector
Dependency ranges changed: mwc-random, primitive, vector, vector-algorithms
Files
- Statistics/Constants.hs +25/−1
- Statistics/Math.hs +173/−42
- statistics.cabal +5/−5
Statistics/Constants.hs view
@@ -15,9 +15,13 @@ , m_huge , m_1_sqrt_2 , m_2_sqrt_pi+ , m_ln_sqrt_2_pi , m_max_exp , m_sqrt_2 , m_sqrt_2_pi+ , m_pos_inf+ , m_neg_inf+ , m_NaN ) where -- | A very large number.@@ -50,7 +54,27 @@ m_1_sqrt_2 = 0.7071067811865475244008443621048490392848359376884740365883 {-# INLINE m_1_sqrt_2 #-} --- | The smallest 'Double' larger than 1.+-- | The smallest 'Double' ε such that 1 + ε ≠ 1. m_epsilon :: Double m_epsilon = encodeFloat (signif+1) expo - 1.0 where (signif,expo) = decodeFloat (1.0::Double)++-- | @log(sqrt((2*pi)) / 2@+m_ln_sqrt_2_pi :: Double+m_ln_sqrt_2_pi = 0.9189385332046727417803297364056176398613974736377834128171+{-# INLINE m_ln_sqrt_2_pi #-}++-- | Positive infinity.+m_pos_inf :: Double+m_pos_inf = 1/0+{-# INLINE m_pos_inf #-}++-- | Negative infinity.+m_neg_inf :: Double+m_neg_inf = -1/0+{-# INLINE m_neg_inf #-}++-- | Not a number.+m_NaN :: Double+m_NaN = 0/0+{-# INLINE m_NaN #-}
Statistics/Math.hs view
@@ -14,57 +14,106 @@ module Statistics.Math ( -- * Functions- chebyshev- , choose- -- ** Factorial functions+ choose+ -- ** Beta function+ , logBeta+ -- ** Chebyshev polynomials+ -- $chebyshev+ , chebyshev+ , chebyshevBroucke+ -- ** Factorial , factorial , logFactorial- -- ** Gamma functions+ -- ** Gamma function , incompleteGamma , logGamma , logGammaL+ -- ** Logarithm+ , log1p -- * References -- $references ) where -import Data.Vector.Generic ((!))+import Data.Int (Int64) import Data.Word (Word64)-import Statistics.Constants (m_sqrt_2_pi)+import Statistics.Constants (m_epsilon, m_sqrt_2_pi, m_ln_sqrt_2_pi, m_NaN,+ m_neg_inf, m_pos_inf) import Statistics.Distribution (cumulative) import Statistics.Distribution.Normal (standard) import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Generic as G +-- $chebyshev+--+-- A Chebyshev polynomial of the first kind is defined by the+-- following recurrence:+--+-- > t 0 _ = 1+-- > t 1 x = x+-- > t n x = 2 * x * t (n-1) x - t (n-2) x+ data C = C {-# UNPACK #-} !Double {-# UNPACK #-} !Double --- | Evaluate a series of Chebyshev polynomials. Uses Clenshaw's--- algorithm.+-- | Evaluate a Chebyshev polynomial of the first kind. Uses+-- Clenshaw's algorithm. chebyshev :: (G.Vector v Double) => Double -- ^ Parameter of each function. -> v Double -- ^ Coefficients of each polynomial term, in increasing order. -> Double-chebyshev x a = fini . U.foldl' step (C 0 0) $ U.enumFromStepN (len - 1) (-1) (len - 1)- where step (C b1 b2) k = C ((a ! k) + x2 * b1 - b2) b1- fini (C b1 b2) = (a ! 0) + x * b1 - b2+chebyshev x a = fini . G.foldr' step (C 0 0) . G.tail $ a+ where step k (C b0 b1) = C (k + x2 * b0 - b1) b0+ fini (C b0 b1) = G.head a + x * b0 - b1 x2 = x * 2- len = G.length a {-# INLINE chebyshev #-} --- | The binomial coefficient.+data B = B {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# UNPACK #-} !Double++-- | Evaluate a Chebyshev polynomial of the first kind. Uses Broucke's+-- ECHEB algorithm, and his convention for coefficient handling, and so+-- gives different results than 'chebyshev' for the same inputs.+chebyshevBroucke :: (G.Vector v Double) =>+ Double -- ^ Parameter of each function.+ -> v Double -- ^ Coefficients of each polynomial term, in increasing order.+ -> Double+chebyshevBroucke x = fini . G.foldr' step (B 0 0 0)+ where step k (B b0 b1 _) = B (k + x2 * b0 - b1) b0 b1+ fini (B b0 _ b2) = (b0 - b2) * 0.5+ x2 = x * 2+{-# INLINE chebyshevBroucke #-}++-- | Quickly compute the natural logarithm of /n/ @`choose`@ /k/, with+-- no checking.+logChooseFast :: Double -> Double -> Double+logChooseFast n k = -log (n + 1) - logBeta (n - k + 1) (k + 1)++-- | Compute the binomial coefficient /n/ @\``choose`\`@ /k/. For+-- values of /k/ > 30, this uses an approximation for performance+-- reasons. The approximation is accurate to 7 decimal places in the+-- worst case, but is typically accurate to 9 decimal places or+-- better. --+-- Example:+-- -- > 7 `choose` 3 == 35 choose :: Int -> Int -> Double n `choose` k- | k > n = 0- | k < 30 = U.foldl' go 1 . U.enumFromTo 1 $ k'- | otherwise = exp $ lg (n+1) - lg (k+1) - lg (n-k+1)- where go a i = a * (nk + j) / j- where j = fromIntegral i :: Double- k' | k > n `div` 2 = n - k- | otherwise = k- nk = fromIntegral (n - k')- lg = logGamma . fromIntegral-{-# INLINE choose #-}+ | k > n = 0+ | k < 30 = U.foldl' go 1 . U.enumFromTo 1 $ k'+ | approx < max64 = fromIntegral . round64 $ approx+ | otherwise = approx+ where+ approx = exp $ logChooseFast (fromIntegral n) (fromIntegral k)+ -- Less numerically stable:+ -- exp $ lg (n+1) - lg (k+1) - lg (n-k+1)+ -- where lg = logGamma . fromIntegral+ go a i = a * (nk + j) / j+ where j = fromIntegral i :: Double+ k' | n_k < k = n_k+ | otherwise = k+ where n_k = n - k+ nk = fromIntegral (n - k')+ max64 = fromIntegral (maxBound :: Int64)+ round64 x = round x :: Int64 data F = F {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 @@ -74,7 +123,7 @@ factorial :: Int -> Double factorial n | n < 0 = error "Statistics.Math.factorial: negative input"- | n <= 1 = 0+ | n <= 1 = 1 | n <= 14 = fini . U.foldl' goLong (F 1 1) $ ns | otherwise = U.foldl' goDouble 1 $ ns where goDouble t k = t * fromIntegral k@@ -82,7 +131,6 @@ where x' = x + 1 fini (F z _) = fromIntegral z ns = U.enumFromTo 2 n-{-# INLINE factorial #-} -- | Compute the natural logarithm of the factorial function. Gives -- 16 decimal digits of precision.@@ -94,18 +142,18 @@ y = 1 / (x * x) z = ((-(5.95238095238e-4 * y) + 7.936500793651e-4) * y - 2.7777777777778e-3) * y + 8.3333333333333e-2-{-# INLINE logFactorial #-} --- | Compute the incomplete gamma integral function γ(/s/,/x/).--- Uses Algorithm AS 239 by Shea.+-- | Compute the normalized lower incomplete gamma function+-- γ(/s/,/x/). Normalization means that+-- γ(∞,/x/)=1. Uses Algorithm AS 239 by Shea. incompleteGamma :: Double -- ^ /s/ -> Double -- ^ /x/ -> Double incompleteGamma x p- | x < 0 || p <= 0 = 1/0+ | x < 0 || p <= 0 = m_pos_inf | x == 0 = 0 | p >= 1000 = norm (3 * sqrt p * ((x/p) ** (1/3) + 1/(9*p) - 1))- | x >= 1e8 = 0+ | x >= 1e8 = 1 | x <= 1 || x < p = let a = p * log x - x - logGamma (p + 1) g = a + log (pearson p 1 1) in if g > limit then exp g else 0@@ -153,7 +201,7 @@ -- ≤ 1e305). logGamma :: Double -> Double logGamma x- | x <= 0 = 1/0+ | x <= 0 = m_pos_inf | x < 1.5 = a + c * ((((r1_4 * b + r1_3) * b + r1_2) * b + r1_1) * b + r1_0) / ((((b + r1_8) * b + r1_7) * b + r1_6) * b + r1_5)@@ -178,19 +226,19 @@ x1 = 1 / x x2 = x1 * x1 - r1_0 = -2.66685511495; r1_1 = -24.4387534237; r1_2 = -21.9698958928- r1_3 = 11.1667541262; r1_4 = 3.13060547623; r1_5 = 0.607771387771- r1_6 = 11.9400905721; r1_7 = 31.4690115749; r1_8 = 15.2346874070+ r1_0 = -2.66685511495; r1_1 = -24.4387534237; r1_2 = -21.9698958928+ r1_3 = 11.1667541262; r1_4 = 3.13060547623; r1_5 = 0.607771387771+ r1_6 = 11.9400905721; r1_7 = 31.4690115749; r1_8 = 15.2346874070 - r2_0 = -78.3359299449; r2_1 = -142.046296688; r2_2 = 137.519416416- r2_3 = 78.6994924154; r2_4 = 4.16438922228; r2_5 = 47.0668766060- r2_6 = 313.399215894; r2_7 = 263.505074721; r2_8 = 43.3400022514+ r2_0 = -78.3359299449; r2_1 = -142.046296688; r2_2 = 137.519416416+ r2_3 = 78.6994924154; r2_4 = 4.16438922228; r2_5 = 47.0668766060+ r2_6 = 313.399215894; r2_7 = 263.505074721; r2_8 = 43.3400022514 - r3_0 = -2.12159572323; r3_1 = 2.30661510616; r3_2 = 2.74647644705- r3_3 = -4.02621119975; r3_4 = -2.29660729780; r3_5 = -1.16328495004- r3_6 = -1.46025937511; r3_7 = -2.42357409629; r3_8 = -5.70691009324+ r3_0 = -2.12159572323e5; r3_1 = 2.30661510616e5; r3_2 = 2.74647644705e4+ r3_3 = -4.02621119975e4; r3_4 = -2.29660729780e3; r3_5 = -1.16328495004e5+ r3_6 = -1.46025937511e5; r3_7 = -2.42357409629e4; r3_8 = -5.70691009324e2 - r4_0 = 0.279195317918525; r4_1 = 0.4917317610505968;+ r4_0 = 0.279195317918525; r4_1 = 0.4917317610505968; r4_2 = 0.0692910599291889; r4_3 = 3.350343815022304 r4_4 = 6.012459259764103 @@ -207,7 +255,7 @@ -- ≤ 1e305). logGammaL :: Double -> Double logGammaL x- | x <= 0 = 1/0+ | x <= 0 = m_pos_inf | otherwise = fini . U.foldl' go (L 0 (x+7)) $ a where fini (L l _) = log (l+a0) + log m_sqrt_2_pi - x65 + (x-0.5) * log x65 go (L l t) k = L (l + k / t) (t-1)@@ -223,7 +271,90 @@ , 676.5203681218835 ] +-- | Compute the log gamma correction factor for @x@ ≥ 10. This+-- correction factor is suitable for an alternate (but less+-- numerically accurate) definition of 'logGamma':+--+-- >lgg x = 0.5 * log(2*pi) + (x-0.5) * log x - x + logGammaCorrection x+logGammaCorrection :: Double -> Double+logGammaCorrection x+ | x < 10 = m_NaN+ | x < big = chebyshevBroucke (t * t * 2 - 1) coeffs / x+ | otherwise = 1 / (x * 12)+ where+ big = 94906265.62425156+ t = 10 / x+ coeffs = U.fromList [+ 0.1666389480451863247205729650822e+0,+ -0.1384948176067563840732986059135e-4,+ 0.9810825646924729426157171547487e-8,+ -0.1809129475572494194263306266719e-10,+ 0.6221098041892605227126015543416e-13,+ -0.3399615005417721944303330599666e-15,+ 0.2683181998482698748957538846666e-17+ ]++-- | Compute the natural logarithm of the beta function.+logBeta :: Double -> Double -> Double+logBeta a b+ | p < 0 = m_NaN+ | p == 0 = m_pos_inf+ | p >= 10 = log q * (-0.5) + m_ln_sqrt_2_pi + logGammaCorrection p + c ++ (p - 0.5) * log ppq + q * log1p(-ppq)+ | q >= 10 = logGamma p + c + p - p * log pq + (q - 0.5) * log1p(-ppq)+ | otherwise = logGamma p + logGamma q - logGamma pq+ where+ p = min a b+ q = max a b+ ppq = p / pq+ pq = p + q+ c = logGammaCorrection q - logGammaCorrection pq++-- | Compute the natural logarithm of 1 + @x@. This is accurate even+-- for values of @x@ near zero, where use of @log(1+x)@ would lose+-- precision.+log1p :: Double -> Double+log1p x+ | x == 0 = 0+ | x == -1 = m_neg_inf+ | x < -1 = m_NaN+ | x' < m_epsilon * 0.5 = x+ | (x >= 0 && x < 1e-8) || (x >= -1e-9 && x < 0)+ = x * (1 - x * 0.5)+ | x' < 0.375 = x * (1 - x * chebyshevBroucke (x / 0.375) coeffs)+ | otherwise = log (1 + x)+ where+ x' = abs x+ coeffs = U.fromList [+ 0.10378693562743769800686267719098e+1,+ -0.13364301504908918098766041553133e+0,+ 0.19408249135520563357926199374750e-1,+ -0.30107551127535777690376537776592e-2,+ 0.48694614797154850090456366509137e-3,+ -0.81054881893175356066809943008622e-4,+ 0.13778847799559524782938251496059e-4,+ -0.23802210894358970251369992914935e-5,+ 0.41640416213865183476391859901989e-6,+ -0.73595828378075994984266837031998e-7,+ 0.13117611876241674949152294345011e-7,+ -0.23546709317742425136696092330175e-8,+ 0.42522773276034997775638052962567e-9,+ -0.77190894134840796826108107493300e-10,+ 0.14075746481359069909215356472191e-10,+ -0.25769072058024680627537078627584e-11,+ 0.47342406666294421849154395005938e-12,+ -0.87249012674742641745301263292675e-13,+ 0.16124614902740551465739833119115e-13,+ -0.29875652015665773006710792416815e-14,+ 0.55480701209082887983041321697279e-15,+ -0.10324619158271569595141333961932e-15+ ]+ -- $references+--+-- * Broucke, R. (1973) Algorithm 446: Ten subroutines for the+-- manipulation of Chebyshev series. /Communications of the ACM/+-- 16(4):254–256. <http://doi.acm.org/10.1145/362003.362037> -- -- * Clenshaw, C.W. (1962) Chebyshev series for mathematical -- functions. /National Physical Laboratory Mathematical Tables 5/,
statistics.cabal view
@@ -1,5 +1,5 @@ name: statistics-version: 0.6.0.0+version: 0.6.0.1 synopsis: A library of statistical types, data, and functions description: This library provides a number of common functions and types useful@@ -54,11 +54,11 @@ build-depends: base < 5, erf,- mwc-random >= 0.5.0.0,- primitive,+ mwc-random >= 0.5.1.4,+ primitive >= 0.3, time,- vector >= 0.5,- vector-algorithms >= 0.3+ vector >= 0.6.0.2,+ vector-algorithms >= 0.3.2 if impl(ghc >= 6.10) build-depends: base >= 4