statistics 0.8.0.5 → 0.9.0.0

raw patch · 14 files changed

+178/−70 lines, 14 filesdep +aesondep +deepseqdep +monad-pardep ~mwc-randomPVP ok

version bump matches the API change (PVP)

Dependencies added: aeson, deepseq, monad-par

Dependency ranges changed: mwc-random

API changes (from Hackage documentation)

Files

• Statistics/Distribution/Binomial.hs +10/−7
• Statistics/Distribution/ChiSquared.hs +6/−2
• Statistics/Distribution/Exponential.hs +6/−3
• Statistics/Distribution/Gamma.hs +9/−2
• Statistics/Distribution/Geometric.hs +5/−2
• Statistics/Distribution/Hypergeometric.hs +9/−8
• Statistics/Distribution/Normal.hs +8/−7
• Statistics/Distribution/Poisson.hs +5/−2
• Statistics/Function.hs +22/−3
• Statistics/KernelDensity.hs +1/−1
• Statistics/Resampling.hs +46/−22
• Statistics/Resampling/Bootstrap.hs +41/−4
• Statistics/Sample.hs +4/−4
• statistics.cabal +6/−3

@@ -23,7 +23,6 @@     , bdProbability     ) where -import Control.Exception (assert) import Data.Typeable (Typeable) import qualified Statistics.Distribution as D import Statistics.Math (choose)@@ -49,7 +48,7 @@     mean = mean  --- This could be slow for bin n+-- This could be slow for big n probability :: BinomialDistribution -> Int -> Double probability (BD n p) k    | k < 0 || k > n = 0@@ -77,12 +76,16 @@ variance (BD n p) = fromIntegral n * p * (1 - p) {-# INLINE variance #-} --- | Construct binomial distribution+-- | Construct binomial distribution. Number of trials must be+--   positive and probability must be in [0,1] range binomial :: Int                 -- ^ Number of trials.          -> Double              -- ^ Probability.          -> BinomialDistribution-binomial n p =-    assert (n > 0) .-    assert (p > 0 && p < 1) $-    BD n p+binomial n p +  | n <= 0         = +    error $ msg ++ "number of trials must be positive. Got " ++ show n+  | p < 0 || p > 1 = +    error $ msg++"probability must be in [0,1] range. Got " ++ show p+  | otherwise      = BD n p+    where msg = "Statistics.Distribution.Binomial.binomial: " {-# INLINE binomial #-}

@@ -34,9 +34,13 @@ chiSquaredNDF (ChiSquared ndf) = ndf {-# INLINE chiSquaredNDF #-} --- | Construct chi-squared distribution. Number of degrees of free+-- | Construct chi-squared distribution. Number of degrees of freedom+--   must be positive. chiSquared :: Int -> ChiSquared-chiSquared x = ChiSquared x+chiSquared n+  | n <= 0    = error $ +     "Statistics.Distribution.ChiSquared.chiSquared: N.D.F. must be positive. Got " ++ show n+  | otherwise = ChiSquared n {-# INLINE chiSquared #-}  instance D.Distribution ChiSquared where

@@ -61,14 +61,17 @@ quantile (ED l) p = -log (1 - p) / l {-# INLINE quantile #-} --- | Create exponential distribution+-- | Create an exponential distribution. exponential :: Double            -- ^ &#955; (scale) parameter.             -> ExponentialDistribution-exponential = ED+exponential l+  | l <= 0 = +    error $ "Statistics.Distribution.Exponential.exponential: scale parameter must be positive. Got " ++ show l+  | otherwise = ED l {-# INLINE exponential #-}  -- | Create exponential distribution from sample. No tests are made to--- check whether it really exponential+-- check whether it truly is exponential. exponentialFromSample :: Sample -> ExponentialDistribution exponentialFromSample = ED . S.mean {-# INLINE exponentialFromSample #-}

@@ -35,8 +35,15 @@     , gdScale :: {-# UNPACK #-} !Double -- ^ Scale parameter, &#977;.     } deriving (Eq, Read, Show, Typeable) -gammaDistr :: Double -> Double -> GammaDistribution-gammaDistr = GD+-- | Create gamma distrivution. Both shape and scale parameters must be positive.+gammaDistr :: Double            -- ^ Shape parameter. /k/+           -> Double            -- ^ Scale parameter, &#977;.+           -> GammaDistribution+gammaDistr k theta+  | k     <= 0 = error $ msg ++ "shape must be positive. Got " ++ show k+  | theta <= 0 = error $ msg ++ "scale must be positive. Got " ++ show theta+  | otherwise  = GD k theta+    where msg = "Statistics.Distribution.Gamma.gammaDistr: " {-# INLINE gammaDistr #-}  instance D.Distribution GammaDistribution where

@@ -51,8 +51,11 @@ -- | Create geometric distribution geometric :: Double                -- ^ Success rate           -> GeometricDistribution-geometric x = assert (x >= 0 && x <= 1)-              GD x+geometric x+  | x < 0 || x > 1 = +    error $ "Statistics.Distribution.Geometric.geometric: probability must be in [0,1] range. Got " ++ show x+  | otherwise      =+    GD x {-# INLINE geometric #-}  probability :: GeometricDistribution -> Int -> Double

@@ -27,9 +27,8 @@     , hdK     ) where -import Control.Exception (assert)-import Data.Typeable (Typeable)-import Statistics.Math (choose)+import Data.Typeable     (Typeable)+import Statistics.Math   (choose) import qualified Statistics.Distribution as D  data HypergeometricDistribution = HD {@@ -66,11 +65,13 @@                -> Int               -- ^ /l/                -> Int               -- ^ /k/                -> HypergeometricDistribution-hypergeometric m l k =-    assert (m >= 0 && m <= l) .-    assert (l > 0) .-    assert (k > 0 && k <= l) $-    HD m l k+hypergeometric m l k+  | not (l > 0)            = error $ msg ++ "l must be positive"+  | not (m >= 0 && m <= l) = error $ msg ++ "m must lie in [0,l] range"+  | not (k > 0 && k <= l)  = error $ msg ++ "k must lie in (0,l] range"+  | otherwise = HD m l k+    where+      msg = "Statistics.Distribution.Hypergeometric.hypergeometric: " {-# INLINE hypergeometric #-}  -- Naive implementation

@@ -20,7 +20,6 @@     , standard     ) where -import Control.Exception (assert) import Data.Number.Erf (erfc) import Data.Typeable (Typeable) import Statistics.Constants (m_sqrt_2, m_sqrt_2_pi)@@ -60,12 +59,14 @@ normalDistr :: Double            -- ^ Mean of distribution             -> Double            -- ^ Variance of distribution             -> NormalDistribution-normalDistr m v = assert (v > 0)-                 ND { mean       = m-                    , variance   = v-                    , ndPdfDenom = m_sqrt_2_pi * sv-                    , ndCdfDenom = m_sqrt_2 * sv-                    }+normalDistr m v+  | v <= 0    = +    error $ "Statistics.Distribution.Normal.normalDistr: variance must be positive. Got " ++ show v+  | otherwise = ND { mean       = m+                   , variance   = v+                   , ndPdfDenom = m_sqrt_2_pi * sv+                   , ndCdfDenom = m_sqrt_2 * sv+                   }     where sv = sqrt v  -- | Create distribution using parameters estimated from

@@ -45,9 +45,12 @@     mean = poissonLambda     {-# INLINE mean #-} --- | Create po+-- | Create poisson distribution. poisson :: Double -> PoissonDistribution-poisson = PD+poisson l+  | l <= 0    = +    error $ "Statistics.Distribution.Poisson.poisson: lambda must be positive. Got " ++ show l+  | otherwise = PD l {-# INLINE poisson #-}  probability :: PoissonDistribution -> Int -> Double

@@ -1,8 +1,7 @@-{-# LANGUAGE Rank2Types #-}-{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE CPP, FlexibleContexts, Rank2Types #-} -- | -- Module    : Statistics.Function--- Copyright : (c) 2009, 2010 Bryan O'Sullivan+-- Copyright : (c) 2009, 2010, 2011 Bryan O'Sullivan -- License   : BSD3 -- -- Maintainer  : bos@serpentine.com@@ -18,12 +17,16 @@     , partialSort     , indexed     , indices+    , nextHighestPowerOfTwo     -- * Vector setup     , create     ) where +#include "MachDeps.h"+ import Control.Exception (assert) import Control.Monad.Primitive (PrimMonad)+import Data.Bits ((.|.), shiftR) import Data.Vector.Generic (modify, unsafeFreeze) import qualified Data.Vector.Algorithms.Intro as I import qualified Data.Vector.Generic as G@@ -74,3 +77,19 @@                                 M.write arr k r                                 loop (k+1) arr {-# INLINE create #-}++-- | Efficiently compute the next highest power of two for a+-- non-negative integer.  If the given value is already a power of+-- two, it is returned unchanged.  If negative, zero is returned.+nextHighestPowerOfTwo :: Int -> Int+nextHighestPowerOfTwo n = o + 1+    where m = n - 1+          o = m+              .|. (m `shiftR` 1)+              .|. (m `shiftR` 2)+              .|. (m `shiftR` 4)+              .|. (m `shiftR` 8)+              .|. (m `shiftR` 16)+#if WORD_SIZE_IN_BITS == 64              +              .|. (m `shiftR` 32)+#endif                

@@ -115,7 +115,7 @@ gaussianKernel :: Kernel gaussianKernel f h p v = exp (-0.5 * u * u) * g     where u = (v - p) / h-          g = f * m_2_sqrt_pi * m_1_sqrt_2+          g = f * 0.5 * m_2_sqrt_pi * m_1_sqrt_2  -- | Kernel density estimator, providing a non-parametric way of -- estimating the PDF of a random variable.

@@ -1,3 +1,5 @@+{-# LANGUAGE BangPatterns #-}+ -- | -- Module    : Statistics.Resampling -- Copyright : (c) 2009, 2010 Bryan O'Sullivan@@ -16,14 +18,16 @@     , resample     ) where -import Control.Monad (forM_, liftM)+import Control.Concurrent (forkIO, newChan, readChan, writeChan)+import Control.Monad (forM_, liftM, replicateM_) import Control.Monad.Primitive (PrimMonad, PrimState) import Data.Vector.Algorithms.Intro (sort) import Data.Vector.Generic (unsafeFreeze)-import Data.Vector.Unboxed ((!))-import Statistics.Function (create, indexed, indices)+import Data.Word (Word32)+import GHC.Conc (numCapabilities)+import Statistics.Function (create, indices) import Statistics.Types (Estimator, Sample)-import System.Random.MWC (Gen, uniform)+import System.Random.MWC (Gen, initialize, uniform, uniformVector) import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Unboxed.Mutable as MU @@ -34,34 +38,54 @@       fromResample :: U.Vector Double     } deriving (Eq, Show) --- | Resample a data set repeatedly, with replacement, computing each--- estimate over the resampled data.-resample :: (PrimMonad m) => Gen (PrimState m) -> [Estimator] -> Int -> Sample -> m [Resample]+-- | /O(e*r*s)/ Resample a data set repeatedly, with replacement,+-- computing each estimate over the resampled data.+--+-- This function is expensive; it has to do work proportional to+-- /e*r*s/, where /e/ is the number of estimation functions, /r/ is+-- the number of resamples to compute, and /s/ is the number of+-- original samples.+--+-- To improve performance, this function will make use of all+-- available CPUs.  At least with GHC 7.0, parallel performance seems+-- best if the parallel garbage collector is disabled (RTS option+-- @-qg@).+resample :: Gen (PrimState IO)+         -> [Estimator]         -- ^ Estimation functions.+         -> Int                 -- ^ Number of resamples to compute.+         -> Sample              -- ^ Original sample.+         -> IO [Resample] resample gen ests numResamples samples = do+  let !numSamples = U.length samples+      ixs = scanl (+) 0 $+            zipWith (+) (replicate numCapabilities q)+                        (replicate r 1 ++ repeat 0)+          where (q,r) = numResamples `quotRem` numCapabilities   results <- mapM (const (MU.new numResamples)) $ ests-  loop 0 (zip ests results)+  done <- newChan+  forM_ (zip ixs (tail ixs)) $ \ (start,!end) -> do+    gen' <- initialize =<< (uniformVector gen 256 :: IO (U.Vector Word32))+    forkIO $ do+      let loop k ers | k >= end = writeChan done ()+                     | otherwise = do+            re <- create numSamples $ \_ -> do+                    r <- uniform gen'+                    return (U.unsafeIndex samples (r `mod` numSamples))+            forM_ ers $ \(est,arr) ->+                MU.write arr k . est $ re+            loop (k+1) ers+      loop start (zip ests results)+  replicateM_ numCapabilities $ readChan done   mapM_ sort results   mapM (liftM Resample . unsafeFreeze) results- where-  loop k ers | k >= numResamples = return ()-             | otherwise = do-    re <- create n $ \_ -> do-            r <- uniform gen-            return (samples ! (abs r `mod` n))-    forM_ ers $ \(est,arr) ->-        MU.write arr k . est $ re-    loop (k+1) ers-  n = U.length samples-{-# INLINE resample #-}  -- | Compute a statistical estimate repeatedly over a sample, each -- time omitting a successive element. jackknife :: Estimator -> Sample -> U.Vector Double jackknife est sample = U.map f . indices $ sample     where f i = est (dropAt i sample)-{-# INLINE jackknife #-}+{- INLINE jackknife #-}  -- | Drop the /k/th element of a vector. dropAt :: U.Unbox e => Int -> U.Vector e -> U.Vector e-dropAt n = U.map snd . U.filter notN . indexed-    where notN (i , _) = i /= n+dropAt n v = U.slice 0 n v U.++ U.slice (n+1) (U.length v - n - 1) v

@@ -1,6 +1,8 @@+{-# LANGUAGE DeriveDataTypeable, OverloadedStrings, RecordWildCards #-}+ -- | -- Module    : Statistics.Resampling.Bootstrap--- Copyright : (c) 2009 Bryan O'Sullivan+-- Copyright : (c) 2009, 2011 Bryan O'Sullivan -- License   : BSD3 -- -- Maintainer  : bos@serpentine.com@@ -13,11 +15,18 @@     (       Estimate(..)     , bootstrapBCA+    , scale     -- * References     -- $references     ) where +import Control.Applicative ((<$>), (<*>), empty)+import Control.DeepSeq (NFData) import Control.Exception (assert)+import Control.Monad.Par (runPar, parMap)+import Data.Aeson.Types+import Data.Data (Data)+import Data.Typeable (Typeable) import Data.Vector.Unboxed ((!)) import Statistics.Distribution (cumulative, quantile) import Statistics.Distribution.Normal@@ -38,8 +47,36 @@     -- the confidence interval).     , estConfidenceLevel :: {-# UNPACK #-} !Double     -- ^ Confidence level of the confidence intervals.-    } deriving (Eq, Show)+    } deriving (Eq, Show, Typeable, Data) +instance NFData Estimate++instance ToJSON Estimate where+    toJSON Estimate{..} = object [+                            "estPoint" .= estPoint+                          , "estLowerBound" .= estLowerBound+                          , "estUpperBound" .= estUpperBound+                          , "estConfidenceLevel" .= estConfidenceLevel+                          ]++instance FromJSON Estimate where+    parseJSON (Object v) = Estimate <$>+                           v .: "estPoint" <*>+                           v .: "estLowerBound" <*>+                           v .: "estUpperBound" <*>+                           v .: "estConfidenceLevel"+    parseJSON _ = empty++-- | Multiply the point, lower bound, and upper bound in an 'Estimate'+-- by the given value.+scale :: Double                 -- ^ Value to multiply by.+      -> Estimate -> Estimate+scale f e@Estimate{..} = e {+                           estPoint = f * estPoint+                         , estLowerBound = f * estLowerBound+                         , estUpperBound = f * estUpperBound+                         }+ estimate :: Double -> Double -> Double -> Double -> Estimate estimate pt lb ub cl =     assert (lb <= ub) .@@ -60,9 +97,9 @@              -> [Estimator]     -- ^ Estimators              -> [Resample]      -- ^ Resampled data              -> [Estimate]-bootstrapBCA confidenceLevel sample =+bootstrapBCA confidenceLevel sample estimators resamples =     assert (confidenceLevel > 0 && confidenceLevel < 1)-    zipWith e+    runPar $ parMap (uncurry e) (zip estimators resamples)   where     e est (Resample resample)       | U.length sample == 1 = estimate pt pt pt confidenceLevel

@@ -60,6 +60,8 @@ -- Operator ^ will be overriden import Prelude hiding ((^)) +-- | Range. The difference between the largest and smallest elements+-- of a sample. range :: (G.Vector v Double) => v Double -> Double range s = hi - lo     where (lo , hi) = minMax s@@ -203,11 +205,9 @@  data V = V {-# UNPACK #-} !Double {-# UNPACK #-} !Double -sqr :: Double -> Double-sqr x = x * x- robustSumVar :: (G.Vector v Double) => Double -> v Double -> Double-robustSumVar m samp = G.sum . G.map (sqr . subtract m) $ samp+robustSumVar m samp = G.sum . G.map (square . subtract m) $ samp+  where square x = x * x {-# INLINE robustSumVar #-}  -- | Maximum likelihood estimate of a sample's variance.  Also known

@@ -1,5 +1,5 @@ name:           statistics-version:        0.8.0.5+version:        0.9.0.0 synopsis:       A library of statistical types, data, and functions description:   This library provides a number of common functions and types useful@@ -26,7 +26,7 @@ bug-reports:    http://bitbucket.org/bos/statistics/issues author:         Bryan O'Sullivan <bos@serpentine.com> maintainer:     Bryan O'Sullivan <bos@serpentine.com>-copyright:      2009, 2010 Bryan O'Sullivan+copyright:      2009, 2010, 2011 Bryan O'Sullivan category:       Math, Statistics build-type:     Simple cabal-version:  >= 1.6@@ -58,9 +58,12 @@   other-modules:     Statistics.Internal   build-depends:+    aeson,     base < 5,+    deepseq >= 1.1.0.2,     erf,-    mwc-random >= 0.8.0.3,+    monad-par >= 0.1.0.1,+    mwc-random >= 0.8.0.5,     primitive >= 0.3,     time,     vector >= 0.7.0.0,