statistics 0.13.1.1 → 0.13.2.0

raw patch · 3 files changed

+72/−7 lines, 3 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• Statistics/Regression.hs +57/−0
• Statistics/Resampling.hs +14/−6
• statistics.cabal +1/−1

@@ -10,15 +10,25 @@       olsRegress     , ols     , rSquare+    , bootstrapRegress     ) where  import Control.Applicative ((<$>))+import Control.Concurrent (forkIO)+import Control.Concurrent.Chan (newChan, readChan, writeChan)+import Control.DeepSeq (rnf)+import Control.Monad (forM_, replicateM)+import GHC.Conc (getNumCapabilities) import Prelude hiding (pred, sum) import Statistics.Function as F import Statistics.Matrix hiding (map) import Statistics.Matrix.Algorithms (qr)+import Statistics.Resampling (splitGen)+import Statistics.Resampling.Bootstrap (Estimate(..)) import Statistics.Sample (mean) import Statistics.Sample.Internal (sum)+import System.Random.MWC (GenIO, uniformR)+import qualified Data.Vector as V import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Unboxed.Mutable as M@@ -97,3 +107,50 @@     r   = sum $ flip U.imap resp $ \i x -> square (x - p i)     t   = sum $ flip U.map resp $ \x -> square (x - mean resp)     p i = sum . flip U.imap coeff $ \j -> (* unsafeIndex pred i j)++-- | Bootstrap a regression function.  Returns both the results of the+-- regression and the requested confidence interval values.+bootstrapRegress :: GenIO+                 -> Int         -- ^ Number of resamples to compute.+                 -> Double      -- ^ Confidence interval.+                 -> ([Vector] -> Vector -> (Vector, Double))+                 -- ^ Regression function.+                 -> [Vector]    -- ^ Predictor vectors.+                 -> Vector      -- ^ Responder vector.+                 -> IO (V.Vector Estimate, Estimate)+bootstrapRegress gen0 numResamples ci rgrss preds0 resp0+  | numResamples < 1   = error $ "bootstrapRegress: number of resamples " +++                                 "must be positive"+  | ci <= 0 || ci >= 1 = error $ "bootstrapRegress: confidence interval " +++                                 "must lie between 0 and 1"+  | otherwise = do+  caps <- getNumCapabilities+  gens <- splitGen caps gen0+  done <- newChan+  forM_ (zip gens (balance caps numResamples)) $ \(gen,count) -> do+    forkIO $ do+      v <- V.replicateM count $ do+           let n = U.length resp0+           ixs <- U.replicateM n $ uniformR (0,n-1) gen+           let resp  = U.backpermute resp0 ixs+               preds = map (flip U.backpermute ixs) preds0+           return $ rgrss preds resp+      rnf v `seq` writeChan done v+  (coeffsv, r2v) <- (G.unzip . V.concat) <$> replicateM caps (readChan done)+  let coeffs  = flip G.imap (G.convert coeffss) $ \i x ->+                est x . U.generate numResamples $ \k -> ((coeffsv G.! k) G.! i)+      r2      = est r2s (G.convert r2v)+      (coeffss, r2s) = rgrss preds0 resp0+      est s v = Estimate s (w G.! lo) (w G.! hi) ci+        where w  = F.sort v+              lo = round c+              hi = truncate (n - c)+              n  = fromIntegral numResamples+              c  = n * ((1 - ci) / 2)+  return (coeffs, r2)++-- | Balance units of work across workers.+balance :: Int -> Int -> [Int]+balance numSlices numItems = zipWith (+) (replicate numSlices q)+                                         (replicate r 1 ++ repeat 0)+ where (q,r) = numItems `quotRem` numSlices

@@ -21,12 +21,12 @@     , jackknifeStdDev     , resample     , estimate+    , splitGen     ) where  import Data.Aeson (FromJSON, ToJSON) import Control.Concurrent (forkIO, newChan, readChan, writeChan)-import Control.Monad (forM_, liftM, replicateM_)-import Control.Monad.Primitive (PrimState)+import Control.Monad (forM_, liftM, replicateM, replicateM_) import Data.Binary (Binary(..)) import Data.Data (Data, Typeable) import Data.Vector.Algorithms.Intro (sort)@@ -39,7 +39,7 @@ import Statistics.Function (indices) import Statistics.Sample (mean, stdDev, variance, varianceUnbiased) import Statistics.Types (Estimator(..), Sample)-import System.Random.MWC (Gen, initialize, uniform, uniformVector)+import System.Random.MWC (GenIO, initialize, uniform, uniformVector) import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Unboxed.Mutable as MU@@ -70,7 +70,7 @@ -- 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)+resample :: GenIO          -> [Estimator]         -- ^ Estimation functions.          -> Int                 -- ^ Number of resamples to compute.          -> Sample              -- ^ Original sample.@@ -83,8 +83,8 @@           where (q,r) = numResamples `quotRem` numCapabilities   results <- mapM (const (MU.new numResamples)) ests   done <- newChan-  forM_ (zip ixs (tail ixs)) $ \ (start,!end) -> do-    gen' <- initialize =<< (uniformVector gen 256 :: IO (U.Vector Word32))+  gens <- splitGen numCapabilities gen+  forM_ (zip3 ixs (tail ixs) gens) $ \ (start,!end,gen') -> do     forkIO $ do       let loop k ers | k >= end = writeChan done ()                      | otherwise = do@@ -175,3 +175,11 @@ singletonErr :: String -> a singletonErr func = error $                     "Statistics.Resampling." ++ func ++ ": singleton input"++-- | Split a generator into several that can run independently.+splitGen :: Int -> GenIO -> IO [GenIO]+splitGen n gen+  | n <= 0    = return []+  | otherwise =+  fmap (gen:) . replicateM (n-1) $+  initialize =<< (uniformVector gen 256 :: IO (U.Vector Word32))

@@ -1,5 +1,5 @@ name:           statistics-version:        0.13.1.1+version:        0.13.2.0 synopsis:       A library of statistical types, data, and functions description:   This library provides a number of common functions and types useful