diff --git a/README b/README --- a/README +++ b/README @@ -8,8 +8,35 @@ estimates for the algorithms used. -Source code +Performance ----------- + +This library has been carefully optimised for high performance. To +obtain the best runtime efficiency, it is imperative to compile +libraries and applications that use this library using a high level of +optimisation. + +Suggested GHC options: + + -O -fvia-C -funbox-strict-fields + +To illustrate, here are the times (in seconds) to generate and sum 250 +million random Word32 values, on a laptop with a 2.4GHz Core2 Duo +P8600 processor, running Fedora 11 and GHC 6.10.3: + + no flags 200+ + -O 1.249 + -O -fvia-C 0.991 + +As the numbers above suggest, compiling without optimisation will +yield unacceptable performance. + + +Get involved! +------------- + +Please feel welcome to contribute new code or bug fixes. You can +fetch the source repository from here: darcs get http://darcs.serpentine.com/statistics diff --git a/Statistics/Distribution.hs b/Statistics/Distribution.hs --- a/Statistics/Distribution.hs +++ b/Statistics/Distribution.hs @@ -21,16 +21,16 @@ -- | The interface shared by all probability distributions. class Distribution d where -- | Probability density function. The probability that a - -- stochastic variable /x/ has the value /X/, i.e. P(/x/=/X/). - probability :: d -> Double -> Double + -- the random variable /X/ has the value /x/, i.e. P(/X/=/x/). + density :: d -> Double -> Double -- | Cumulative distribution function. The probability that a - -- stochastic variable /x/ is less than /X/, i.e. P(/x/ Double -> Double + -- random variable /X/ is less than /x/, i.e. P(/X/≤/x/). + cumulative :: d -> Double -> Double -- | Inverse of the cumulative distribution function. The value - -- /X/ for which P(/x/ Double -> Double + -- /x/ for which P(/X/≤/x/). + quantile :: d -> Double -> Double class Distribution d => Mean d where mean :: d -> Double @@ -38,6 +38,8 @@ class Mean d => Variance d where variance :: d -> Double +data P = P {-# UNPACK #-} !Double {-# UNPACK #-} !Double + -- | Approximate the value of /X/ for which P(/x/>/X/)=/p/. -- -- This method uses a combination of Newton-Raphson iteration and @@ -57,13 +59,13 @@ | otherwise = loop (i+1) dx'' x'' lo' hi' where err = cumulative d x - prob - (lo',hi') | err < 0 = (x, hi) - | otherwise = (lo, x) - pdf = probability d x - (dx',x') | pdf /= 0 = (err / pdf, x - dx) - | otherwise = (dx, x) - (dx'',x'') - | x' < lo' || x' > hi' || pdf == 0 = (x'-x, (lo + hi) / 2) - | otherwise = (dx', x') + P lo' hi' | err < 0 = P x hi + | otherwise = P lo x + pdf = density d x + P dx' x' | pdf /= 0 = P (err / pdf) (x - dx) + | otherwise = P dx x + P dx'' x'' + | x' < lo' || x' > hi' || pdf == 0 = P (x'-x) ((lo + hi) / 2) + | otherwise = P dx' x' accuracy = 1e-15 maxIters = 150 diff --git a/Statistics/Distribution/Binomial.hs b/Statistics/Distribution/Binomial.hs --- a/Statistics/Distribution/Binomial.hs +++ b/Statistics/Distribution/Binomial.hs @@ -38,9 +38,9 @@ } deriving (Eq, Read, Show, Typeable) instance D.Distribution BinomialDistribution where - probability = probability + density = density cumulative = cumulative - inverse = inverse + quantile = quantile instance D.Variance BinomialDistribution where variance = variance @@ -48,16 +48,16 @@ instance D.Mean BinomialDistribution where mean = mean -probability :: BinomialDistribution -> Double -> Double -probability (BD n p) x = +density :: BinomialDistribution -> Double -> Double +density (BD n p) x = fromIntegral (n `choose` floor x) * p ** x * (1-p) ** (fromIntegral n-x) cumulative :: BinomialDistribution -> Double -> Double cumulative d = - sumU . mapU (probability d . fromIntegral) . enumFromToU (0::Int) . floor + sumU . mapU (density d . fromIntegral) . enumFromToU (0::Int) . floor -inverse :: BinomialDistribution -> Double -> Double -inverse d@(BD n _p) p = D.findRoot d p (n'/2) 0 n' +quantile :: BinomialDistribution -> Double -> Double +quantile d@(BD n _p) p = D.findRoot d p (n'/2) 0 n' where n' = fromIntegral n mean :: BinomialDistribution -> Double diff --git a/Statistics/Distribution/Exponential.hs b/Statistics/Distribution/Exponential.hs --- a/Statistics/Distribution/Exponential.hs +++ b/Statistics/Distribution/Exponential.hs @@ -33,12 +33,12 @@ } deriving (Eq, Read, Show, Typeable) instance D.Distribution ExponentialDistribution where - probability (ED l) x = l * exp (-l * x) - {-# INLINE probability #-} - cumulative (ED l) x = 1 - exp (-l * x) + density (ED l) x = l * exp (-l * x) + {-# INLINE density #-} + cumulative (ED l) x = 1 - exp (-l * x) {-# INLINE cumulative #-} - inverse (ED l) p = -log (1 - p) / l - {-# INLINE inverse #-} + quantile (ED l) p = -log (1 - p) / l + {-# INLINE quantile #-} instance D.Variance ExponentialDistribution where variance (ED l) = 1 / (l * l) diff --git a/Statistics/Distribution/Gamma.hs b/Statistics/Distribution/Gamma.hs --- a/Statistics/Distribution/Gamma.hs +++ b/Statistics/Distribution/Gamma.hs @@ -38,9 +38,9 @@ } deriving (Eq, Read, Show, Typeable) instance D.Distribution GammaDistribution where - probability = probability - cumulative = cumulative - inverse = inverse + density = density + cumulative = cumulative + quantile = quantile instance D.Variance GammaDistribution where variance (GD a l) = a / (l * l) @@ -50,17 +50,17 @@ mean (GD a l) = a / l {-# INLINE mean #-} -probability :: GammaDistribution -> Double -> Double -probability (GD a l) x = x ** (a-1) * exp (-x/l) / (exp (logGamma a) * l ** a) -{-# INLINE probability #-} +density :: GammaDistribution -> Double -> Double +density (GD a l) x = x ** (a-1) * exp (-x/l) / (exp (logGamma a) * l ** a) +{-# INLINE density #-} cumulative :: GammaDistribution -> Double -> Double cumulative (GD a l) x = incompleteGamma a (x/l) / exp (logGamma a) {-# INLINE cumulative #-} -inverse :: GammaDistribution -> Double -> Double -inverse d p +quantile :: GammaDistribution -> Double -> Double +quantile d p | p == 0 = -1/0 | p == 1 = 1/0 | otherwise = D.findRoot d p (gdShape d) 0 m_huge -{-# INLINE inverse #-} +{-# INLINE quantile #-} diff --git a/Statistics/Distribution/Geometric.hs b/Statistics/Distribution/Geometric.hs --- a/Statistics/Distribution/Geometric.hs +++ b/Statistics/Distribution/Geometric.hs @@ -35,9 +35,9 @@ } deriving (Eq, Read, Show, Typeable) instance D.Distribution GeometricDistribution where - probability = probability - cumulative = cumulative - inverse = inverse + density = density + cumulative = cumulative + quantile = quantile instance D.Variance GeometricDistribution where variance (GD s) = (1 - s) / (s * s) @@ -52,14 +52,14 @@ GD x {-# INLINE fromSuccess #-} -probability :: GeometricDistribution -> Double -> Double -probability (GD s) x = s * (1-s) ** (x-1) -{-# INLINE probability #-} +density :: GeometricDistribution -> Double -> Double +density (GD s) x = s * (1-s) ** (x-1) +{-# INLINE density #-} cumulative :: GeometricDistribution -> Double -> Double cumulative (GD s) x = 1 - (1-s) ** x {-# INLINE cumulative #-} -inverse :: GeometricDistribution -> Double -> Double -inverse (GD s) p = log (1 - p) / log (1 - s) -{-# INLINE inverse #-} +quantile :: GeometricDistribution -> Double -> Double +quantile (GD s) p = log (1 - p) / log (1 - s) +{-# INLINE quantile #-} diff --git a/Statistics/Distribution/Hypergeometric.hs b/Statistics/Distribution/Hypergeometric.hs --- a/Statistics/Distribution/Hypergeometric.hs +++ b/Statistics/Distribution/Hypergeometric.hs @@ -41,9 +41,9 @@ } deriving (Eq, Read, Show, Typeable) instance D.Distribution HypergeometricDistribution where - probability = probability - cumulative = cumulative - inverse = inverse + density = density + cumulative = cumulative + quantile = quantile instance D.Variance HypergeometricDistribution where variance = variance @@ -74,8 +74,8 @@ HD m l k {-# INLINE fromParams #-} -probability :: HypergeometricDistribution -> Double -> Double -probability (HD mi li ki) x +density :: HypergeometricDistribution -> Double -> Double +density (HD mi li ki) x | l <= 70 = (mi <> xi) * ((li - mi) <> (ki - xi)) / (li <> ki) | r > maxVal = 1/0 | otherwise = exp r @@ -89,7 +89,7 @@ m = fromIntegral mi l = fromIntegral li k = fromIntegral ki -{-# INLINE probability #-} +{-# INLINE density #-} cumulative :: HypergeometricDistribution -> Double -> Double cumulative d@(HD m l k) x @@ -99,9 +99,9 @@ where imin = max 0 (k - l + m) imax = min k m - r = sumU . mapU (probability d . fromIntegral) . enumFromToU imin . floor $ x + r = sumU . mapU (density d . fromIntegral) . enumFromToU imin . floor $ x {-# INLINE cumulative #-} -inverse :: HypergeometricDistribution -> Double -> Double -inverse = error "Statistics.Distribution.Hypergeometric.inverse: not yet implemented" -{-# INLINE inverse #-} +quantile :: HypergeometricDistribution -> Double -> Double +quantile = error "Statistics.Distribution.Hypergeometric.quantile: not yet implemented" +{-# INLINE quantile #-} diff --git a/Statistics/Distribution/Normal.hs b/Statistics/Distribution/Normal.hs --- a/Statistics/Distribution/Normal.hs +++ b/Statistics/Distribution/Normal.hs @@ -37,9 +37,9 @@ } deriving (Eq, Read, Show, Typeable) instance D.Distribution NormalDistribution where - probability = probability - cumulative = cumulative - inverse = inverse + density = density + cumulative = cumulative + quantile = quantile instance D.Variance NormalDistribution where variance = variance @@ -68,15 +68,15 @@ fromSample :: Sample -> NormalDistribution fromSample a = fromParams (S.mean a) (S.variance a) -probability :: NormalDistribution -> Double -> Double -probability d x = exp (-xm * xm / (2 * variance d)) / ndPdfDenom d +density :: NormalDistribution -> Double -> Double +density d x = exp (-xm * xm / (2 * variance d)) / ndPdfDenom d where xm = x - mean d cumulative :: NormalDistribution -> Double -> Double cumulative d x = erfc (-(x-mean d) / ndCdfDenom d) / 2 -inverse :: NormalDistribution -> Double -> Double -inverse d p +quantile :: NormalDistribution -> Double -> Double +quantile d p | p == 0 = -m_huge | p == 1 = m_huge | p == 0.5 = mean d diff --git a/Statistics/Distribution/Poisson.hs b/Statistics/Distribution/Poisson.hs --- a/Statistics/Distribution/Poisson.hs +++ b/Statistics/Distribution/Poisson.hs @@ -32,9 +32,9 @@ } deriving (Eq, Read, Show, Typeable) instance D.Distribution PoissonDistribution where - probability = probability - cumulative = cumulative - inverse = inverse + density = density + cumulative = cumulative + quantile = quantile instance D.Variance PoissonDistribution where variance = pdLambda @@ -48,16 +48,16 @@ fromLambda = PD {-# INLINE fromLambda #-} -probability :: PoissonDistribution -> Double -> Double -probability (PD l) x = exp (x * log l - l - logGamma x) -{-# INLINE probability #-} +density :: PoissonDistribution -> Double -> Double +density (PD l) x = exp (x * log l - l - logGamma x) +{-# INLINE density #-} cumulative :: PoissonDistribution -> Double -> Double -cumulative d = sumU . mapU (probability d . fromIntegral) . +cumulative d = sumU . mapU (density d . fromIntegral) . enumFromToU (0::Int) . floor {-# INLINE cumulative #-} -inverse :: PoissonDistribution -> Double -> Double -inverse d p = fromIntegral . r $ D.findRoot d p (pdLambda d) 0 m_huge +quantile :: PoissonDistribution -> Double -> Double +quantile d p = fromIntegral . r $ D.findRoot d p (pdLambda d) 0 m_huge where r = round :: Double -> Int -{-# INLINE inverse #-} +{-# INLINE quantile #-} diff --git a/Statistics/Function.hs b/Statistics/Function.hs --- a/Statistics/Function.hs +++ b/Statistics/Function.hs @@ -1,4 +1,4 @@ -{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE Rank2Types, TypeOperators #-} -- | -- Module : Statistics.Function -- Copyright : (c) 2009 Bryan O'Sullivan @@ -19,7 +19,7 @@ ) where import Control.Exception (assert) -import Control.Monad.ST (unsafeSTToIO) +import Control.Monad.ST (ST) import Data.Array.Vector.Algorithms.Combinators (apply) import Data.Array.Vector import qualified Data.Array.Vector.Algorithms.Intro as I @@ -49,12 +49,12 @@ {-# INLINE minMax #-} -- | Create an array, using the given action to populate each element. -createU :: (UA e) => Int -> (Int -> IO e) -> IO (UArr e) +createU :: (UA e) => forall s. Int -> (Int -> ST s e) -> ST s (UArr e) createU size itemAt = assert (size >= 0) $ - unsafeSTToIO (newMU size) >>= loop 0 + newMU size >>= loop 0 where - loop k arr | k >= size = unsafeSTToIO (unsafeFreezeAllMU arr) + loop k arr | k >= size = unsafeFreezeAllMU arr | otherwise = do r <- itemAt k - unsafeSTToIO (writeMU arr k r) + writeMU arr k r loop (k+1) arr diff --git a/Statistics/RandomVariate.hs b/Statistics/RandomVariate.hs new file mode 100644 --- /dev/null +++ b/Statistics/RandomVariate.hs @@ -0,0 +1,445 @@ +{-# LANGUAGE BangPatterns, CPP, MagicHash, Rank2Types, ScopedTypeVariables #-} +-- | +-- Module : Statistics.RandomVariate +-- Copyright : (c) 2009 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Pseudo-random variate generation. + +module Statistics.RandomVariate + ( + -- * Types + Gen + , Variate(..) + -- * Other distributions + , normal + -- * Creation + , create + , initialize + , withSystemRandom + -- * Helper functions + , uniformArray + -- * References + -- $references + ) where + +#if defined(__GLASGOW_HASKELL__) && !defined(__HADDOCK__) +#include "MachDeps.h" +#endif + +import Control.Exception (IOException, catch) +import Control.Monad (ap, unless) +import Control.Monad.ST (ST, runST) +import Data.Array.Vector +import Data.Bits ((.&.), (.|.), xor) +import Data.IORef (atomicModifyIORef, newIORef) +import Data.Int (Int8, Int16, Int32, Int64) +import Data.Ratio ((%), numerator) +import Data.Time.Clock.POSIX (getPOSIXTime) +import Data.Word (Word, Word8, Word16, Word32, Word64) +import Foreign.Marshal.Alloc (allocaBytes) +import Foreign.Marshal.Array (peekArray) +import GHC.Base (Int(I#)) +import GHC.Word (Word64(W64#), uncheckedShiftL64#, uncheckedShiftRL64#) +import Prelude hiding (catch) +import System.CPUTime (cpuTimePrecision, getCPUTime) +import System.IO (IOMode(..), hGetBuf, hPutStrLn, stderr, withBinaryFile) +import System.IO.Unsafe (unsafePerformIO) + +-- | The class of types for which we can generate uniformly +-- distributed random variates. +-- +-- The uniform PRNG uses Marsaglia's MWC256 (also known as MWC8222) +-- multiply-with-carry generator, which has a period of 2^8222 and +-- fares well in tests of randomness. It is also extremely fast, +-- between 2 and 3 times faster than the Mersenne Twister. +-- +-- /Note/: Marsaglia's PRNG is not known to be cryptographically +-- secure, so you should not use it for cryptographic operations. +class Variate a where + -- | Generate a single uniformly distributed random variate. The + -- range of values produced varies by type: + -- + -- * For fixed-width integral types, the type's entire range is + -- used. + -- + -- * For floating point numbers, the range (0,1] is used. Zero is + -- explicitly excluded, to allow variates to be used in + -- statistical calculations that require non-zero values + -- (e.g. uses of the 'log' function). + -- + -- * The range of random 'Integer' variates is the same as for + -- 'Int'. + uniform :: Gen s -> ST s a + +-- Thanks to Duncan Coutts for finding the pattern below for +-- strong-arming GHC 6.10's inliner into behaving itself. This makes +-- a 2x difference to performance compared to the following: +-- +-- > uniform = uniform1 fromIntegral + +instance Variate Int8 where + uniform = f where f = uniform1 fromIntegral + {-# INLINE f #-} + +instance Variate Int16 where + uniform = f where f = uniform1 fromIntegral + {-# INLINE f #-} + +instance Variate Int32 where + uniform = f where f = uniform1 fromIntegral + {-# INLINE f #-} + +instance Variate Int64 where + uniform = f where f = uniform2 wordsTo64Bit + {-# INLINE f #-} + +instance Variate Word8 where + uniform = f where f = uniform1 fromIntegral + {-# INLINE f #-} + +instance Variate Word16 where + uniform = f where f = uniform1 fromIntegral + {-# INLINE f #-} + +instance Variate Word32 where + uniform = uniformWord32 + +instance Variate Word64 where + uniform = f where f = uniform2 wordsTo64Bit + {-# INLINE f #-} + +instance Variate Bool where + uniform = f where f = uniform1 wordToBool + {-# INLINE f #-} + +instance Variate Float where + uniform = f where f = uniform1 wordToFloat + {-# INLINE f #-} + +instance Variate Double where + uniform = f where f = uniform2 wordsToDouble + {-# INLINE f #-} + +instance Variate Int where +#if WORD_SIZE_IN_BITS < 64 + uniform = f where f = uniform1 fromIntegral +#else + uniform = f where f = uniform2 wordsTo64Bit +#endif + {-# INLINE f #-} + +instance Variate Word where +#if WORD_SIZE_IN_BITS < 64 + uniform = f where f = uniform1 fromIntegral +#else + uniform = f where f = uniform2 wordsTo64Bit +#endif + {-# INLINE f #-} + +instance Variate Integer where + uniform = f where f g = do + u <- uniform g + return $! fromIntegral (u :: Int) + {-# INLINE f #-} + +instance (Variate a, Variate b) => Variate (a,b) where + uniform = f where f g = (,) `fmap` uniform g `ap` uniform g + {-# INLINE f #-} + +instance (Variate a, Variate b, Variate c) => Variate (a,b,c) where + uniform = f where f g = (,,) `fmap` uniform g `ap` uniform g `ap` uniform g + {-# INLINE f #-} + +instance (Variate a, Variate b, Variate c, Variate d) => Variate (a,b,c,d) where + uniform = f + where f g = (,,,) `fmap` uniform g `ap` uniform g `ap` uniform g + `ap` uniform g + {-# INLINE f #-} + +wordsTo64Bit :: Integral a => Word32 -> Word32 -> a +wordsTo64Bit a b = + fromIntegral ((fromIntegral a `shiftL` 32) .|. fromIntegral b) +{-# INLINE wordsTo64Bit #-} + +wordToBool :: Word32 -> Bool +wordToBool i = (i .&. 1) /= 0 +{-# INLINE wordToBool #-} + +wordToFloat :: Word32 -> Float +wordToFloat x = (fromIntegral i * m_inv_32) + 0.5 + m_inv_33 + where m_inv_33 = 1.16415321826934814453125e-10 + m_inv_32 = 2.3283064365386962890625e-10 + i = fromIntegral x :: Int32 +{-# INLINE wordToFloat #-} + +wordsToDouble :: Word32 -> Word32 -> Double +wordsToDouble x y = (fromIntegral a * m_inv_32 + (0.5 + m_inv_53) + + fromIntegral (b .&. 0xFFFFF) * m_inv_52) + where m_inv_52 = 2.220446049250313080847263336181640625e-16 + m_inv_53 = 1.1102230246251565404236316680908203125e-16 + m_inv_32 = 2.3283064365386962890625e-10 + a = fromIntegral x :: Int32 + b = fromIntegral y :: Int32 +{-# INLINE wordsToDouble #-} + +-- | State of the pseudo-random number generator. +newtype Gen s = Gen (MUArr Word32 s) + +ioff, coff :: Int +ioff = 256 +coff = 257 + +-- | Create a generator for variates using a fixed seed. +create :: ST s (Gen s) +create = initialize defaultSeed +{-# INLINE create #-} + +-- | Create a generator for variates using the given seed, of which up +-- to 256 elements will be used. For arrays of less than 256 +-- elements, part of the default seed will be used to finish +-- initializing the generator's state. +-- +-- Examples: +-- +-- > initialize (singletonU 42) +-- +-- > initialize (toU [4, 8, 15, 16, 23, 42]) +-- +-- If a seed contains fewer than 256 elements, it is first used +-- verbatim, then its elements are 'xor'ed against elements of the +-- default seed until 256 elements are reached. +initialize :: UArr Word32 -> ST s (Gen s) +initialize seed = do + q <- newMU 258 + fill q + writeMU q ioff 255 + writeMU q coff 362436 + return (Gen q) + where fill q = go 0 where + go i | i == 256 = return () + | otherwise = writeMU q i s >> go (i+1) + where s | i >= fini = if fini == 0 + then indexU defaultSeed i + else indexU defaultSeed i `xor` + indexU seed (i `mod` fini) + | otherwise = indexU seed i + fini = lengthU seed +{-# INLINE initialize #-} + +-- | Using the current time as a seed, perform an action that uses a +-- random variate generator. This is a horrible fallback for Windows +-- systems. +withTime :: (forall s. Gen s -> ST s a) -> IO a +withTime act = do + c <- (numerator . (%cpuTimePrecision)) `fmap` getCPUTime + t <- toRational `fmap` getPOSIXTime + let n = fromIntegral (numerator t) :: Word64 + seed = [fromIntegral c, fromIntegral n, fromIntegral (n `shiftR` 32)] + return . runST $ initialize (toU seed) >>= act + +-- | Seed a PRNG with data from the system's fast source of +-- pseudo-random numbers (\"\/dev\/urandom\" on Unix-like systems), +-- then run the given action. +-- +-- /Note/: on Windows, this code does not yet use the native +-- Cryptographic API as a source of random numbers (it uses the system +-- clock instead). As a result, the sequences it generates may not be +-- highly independent. +withSystemRandom :: (forall s. Gen s -> ST s a) -> IO a +withSystemRandom act = tryRandom `catch` \(_::IOException) -> do + seen <- atomicModifyIORef warned ((,) True) + unless seen $ do + hPutStrLn stderr ("Warning: Couldn't open " ++ show random) + hPutStrLn stderr ("Warning: using system clock for seed instead " ++ + "(quality will be lower)") + withTime act + where tryRandom = do + let nbytes = 1024 + ws <- allocaBytes nbytes $ \buf -> do + nread <- withBinaryFile random ReadMode $ + \h -> hGetBuf h buf nbytes + peekArray (nread `div` 4) buf + return . runST $ initialize (toU ws) >>= act + random = "/dev/urandom" + warned = unsafePerformIO $ newIORef False + {-# NOINLINE warned #-} + +-- | Unchecked 64-bit left shift. +shiftL :: Word64 -> Int -> Word64 +shiftL (W64# x#) (I# i#) = W64# (x# `uncheckedShiftL64#` i#) + +-- | Unchecked 64-bit right shift. +shiftR :: Word64 -> Int -> Word64 +shiftR (W64# x#) (I# i#) = W64# (x# `uncheckedShiftRL64#` i#) + +-- | Compute the next index into the state pool. This is simply +-- addition modulo 256. +nextIndex :: Integral a => a -> Int +nextIndex i = fromIntegral j + where j = fromIntegral (i+1) :: Word8 + +uniformWord32 :: Gen s -> ST s Word32 +uniformWord32 (Gen q) = do + let a = 809430660 :: Word64 + i <- nextIndex `fmap` readMU q ioff + c <- fromIntegral `fmap` readMU q coff + qi <- fromIntegral `fmap` readMU q i + let t = a * qi + c + t32 = fromIntegral t + writeMU q i t32 + writeMU q ioff (fromIntegral i) + writeMU q coff (fromIntegral (t `shiftR` 32)) + return t32 +{-# INLINE uniformWord32 #-} + +uniform1 :: (Word32 -> a) -> Gen s -> ST s a +uniform1 f gen = do + i <- uniformWord32 gen + return $! f i +{-# INLINE uniform1 #-} + +uniform2 :: (Word32 -> Word32 -> a) -> Gen s -> ST s a +uniform2 f (Gen q) = do + let a = 809430660 :: Word64 + i <- nextIndex `fmap` readMU q ioff + let j = nextIndex i + c <- fromIntegral `fmap` readMU q coff + qi <- fromIntegral `fmap` readMU q i + qj <- fromIntegral `fmap` readMU q j + let t = a * qi + c + t32 = fromIntegral t + c' = t `shiftR` 32 + u = a * qj + c' + u32 = fromIntegral u + writeMU q i t32 + writeMU q j u32 + writeMU q ioff (fromIntegral j) + writeMU q coff (fromIntegral (u `shiftR` 32)) + return $! f t32 u32 +{-# INLINE uniform2 #-} + +-- | Generate an array of pseudo-random variates. This is not +-- necessarily faster than invoking 'uniform' repeatedly in a loop, +-- but it may be more convenient to use in some situations. +uniformArray :: (UA a, Variate a) => Gen s -> Int -> ST s (UArr a) +uniformArray gen n = newMU n >>= loop + where + loop mu = go 0 + where go !i | i >= n = unsafeFreezeAllMU mu + | otherwise = uniform gen >>= writeMU mu i >> go (i+1) +{-# INLINE uniformArray #-} + +-- | Generate a normally distributed random variate. +-- +-- The implementation uses Doornik's modified ziggurat algorithm. +-- Compared to the ziggurat algorithm usually used, this is slower, +-- but generates more independent variates that pass stringent tests +-- of randomness. +normal :: Gen s -> ST s Double +normal gen = loop + where + loop = do + u <- (subtract 1 . (*2)) `fmap` uniform gen + ri <- uniform gen + let i = fromIntegral ((ri :: Word32) .&. 127) + bi = indexU blocks i + bj = indexU blocks (i+1) + if abs u < indexU ratios i + then return $! u * bi + else if i == 0 + then normalTail (u < 0) + else do + let x = u * bi + xx = x * x + d = exp (-0.5 * (bi * bi - xx)) + e = exp (-0.5 * (bj * bj - xx)) + c <- uniform gen + if e + c * (d - e) < 1 + then return x + else loop + blocks = let f = exp (-0.5 * r * r) + in (`snocU` 0) . consU (v/f) . consU r . unfoldU 126 go $ (r :*: f) + where + go (b :*: g) = JustS (h :*: (h :*: exp (-0.5 * h * h))) + where h = sqrt (-2 * log (v / b + g)) + v = 9.91256303526217e-3 + r = 3.442619855899 + ratios = zipWithU (/) (tailU blocks) blocks + normalTail neg = tailing + where tailing = do + x <- ((/r) . log) `fmap` uniform gen + y <- log `fmap` uniform gen + if y * (-2) < x * x + then tailing + else return $! if neg then x - r else r - x + +defaultSeed :: UArr Word32 +defaultSeed = toU [ + 0x7042e8b3, 0x06f7f4c5, 0x789ea382, 0x6fb15ad8, 0x54f7a879, 0x0474b184, + 0xb3f8f692, 0x4114ea35, 0xb6af0230, 0xebb457d2, 0x47693630, 0x15bc0433, + 0x2e1e5b18, 0xbe91129c, 0xcc0815a0, 0xb1260436, 0xd6f605b1, 0xeaadd777, + 0x8f59f791, 0xe7149ed9, 0x72d49dd5, 0xd68d9ded, 0xe2a13153, 0x67648eab, + 0x48d6a1a1, 0xa69ab6d7, 0x236f34ec, 0x4e717a21, 0x9d07553d, 0x6683a701, + 0x19004315, 0x7b6429c5, 0x84964f99, 0x982eb292, 0x3a8be83e, 0xc1df1845, + 0x3cf7b527, 0xb66a7d3f, 0xf93f6838, 0x736b1c85, 0x5f0825c1, 0x37e9904b, + 0x724cd7b3, 0xfdcb7a46, 0xfdd39f52, 0x715506d5, 0xbd1b6637, 0xadabc0c0, + 0x219037fc, 0x9d71b317, 0x3bec717b, 0xd4501d20, 0xd95ea1c9, 0xbe717202, + 0xa254bd61, 0xd78a6c5b, 0x043a5b16, 0x0f447a25, 0xf4862a00, 0x48a48b75, + 0x1e580143, 0xd5b6a11b, 0x6fb5b0a4, 0x5aaf27f9, 0x668bcd0e, 0x3fdf18fd, + 0x8fdcec4a, 0x5255ce87, 0xa1b24dbf, 0x3ee4c2e1, 0x9087eea2, 0xa4131b26, + 0x694531a5, 0xa143d867, 0xd9f77c03, 0xf0085918, 0x1e85071c, 0x164d1aba, + 0xe61abab5, 0xb8b0c124, 0x84899697, 0xea022359, 0x0cc7fa0c, 0xd6499adf, + 0x746da638, 0xd9e5d200, 0xefb3360b, 0x9426716a, 0xabddf8c2, 0xdd1ed9e4, + 0x17e1d567, 0xa9a65000, 0x2f37dbc5, 0x9a4b8fd5, 0xaeb22492, 0x0ebe8845, + 0xd89dd090, 0xcfbb88c6, 0xb1325561, 0x6d811d90, 0x03aa86f4, 0xbddba397, + 0x0986b9ed, 0x6f4cfc69, 0xc02b43bc, 0xee916274, 0xde7d9659, 0x7d3afd93, + 0xf52a7095, 0xf21a009c, 0xfd3f795e, 0x98cef25b, 0x6cb3af61, 0x6fa0e310, + 0x0196d036, 0xbc198bca, 0x15b0412d, 0xde454349, 0x5719472b, 0x8244ebce, + 0xee61afc6, 0xa60c9cb5, 0x1f4d1fd0, 0xe4fb3059, 0xab9ec0f9, 0x8d8b0255, + 0x4e7430bf, 0x3a22aa6b, 0x27de22d3, 0x60c4b6e6, 0x0cf61eb3, 0x469a87df, + 0xa4da1388, 0xf650f6aa, 0x3db87d68, 0xcdb6964c, 0xb2649b6c, 0x6a880fa9, + 0x1b0c845b, 0xe0af2f28, 0xfc1d5da9, 0xf64878a6, 0x667ca525, 0x2114b1ce, + 0x2d119ae3, 0x8d29d3bf, 0x1a1b4922, 0x3132980e, 0xd59e4385, 0x4dbd49b8, + 0x2de0bb05, 0xd6c96598, 0xb4c527c3, 0xb5562afc, 0x61eeb602, 0x05aa192a, + 0x7d127e77, 0xc719222d, 0xde7cf8db, 0x2de439b8, 0x250b5f1a, 0xd7b21053, + 0xef6c14a1, 0x2041f80f, 0xc287332e, 0xbb1dbfd3, 0x783bb979, 0x9a2e6327, + 0x6eb03027, 0x0225fa2f, 0xa319bc89, 0x864112d4, 0xfe990445, 0xe5e2e07c, + 0xf7c6acb8, 0x1bc92142, 0x12e9b40e, 0x2979282d, 0x05278e70, 0xe160ba4c, + 0xc1de0909, 0x458b9bf4, 0xbfce9c94, 0xa276f72a, 0x8441597d, 0x67adc2da, + 0x6162b854, 0x7f9b2f4a, 0x0d995b6b, 0x193b643d, 0x399362b3, 0x8b653a4b, + 0x1028d2db, 0x2b3df842, 0x6eecafaf, 0x261667e9, 0x9c7e8cda, 0x46063eab, + 0x7ce7a3a1, 0xadc899c9, 0x017291c4, 0x528d1a93, 0x9a1ee498, 0xbb7d4d43, + 0x7837f0ed, 0x34a230cc, 0x614a628d, 0xb03f93b8, 0xd72e3b08, 0x604c98db, + 0x3cfacb79, 0x8b81646a, 0xc0f082fa, 0xd1f92388, 0xe5a91e39, 0xf95c756d, + 0x1177742f, 0xf8819323, 0x5c060b80, 0x96c1cd8f, 0x47d7b440, 0xbbb84197, + 0x35f749cc, 0x95b0e132, 0x8d90ad54, 0x5c3f9423, 0x4994005b, 0xb58f53b9, + 0x32df7348, 0x60f61c29, 0x9eae2f32, 0x85a3d398, 0x3b995dd4, 0x94c5e460, + 0x8e54b9f3, 0x87bc6e2a, 0x90bbf1ea, 0x55d44719, 0x2cbbfe6e, 0x439d82f0, + 0x4eb3782d, 0xc3f1e669, 0x61ff8d9e, 0x0909238d, 0xef406165, 0x09c1d762, + 0x705d184f, 0x188f2cc4, 0x9c5aa12a, 0xc7a5d70e, 0xbc78cb1b, 0x1d26ae62, + 0x23f96ae3, 0xd456bf32, 0xe4654f55, 0x31462bd8 ] + +-- $references +-- +-- * Doornik, J.A. (2005) An improved ziggurat method to generate +-- normal random samples. Mimeo, Nuffield College, University of +-- Oxford. +-- +-- * Doornik, J.A. (2007) Conversion of high-period random numbers to +-- floating point. +-- /ACM Transactions on Modeling and Computer Simulation/ 17(1). +-- +-- +-- * Marsaglia, G. (2003) Seeds for random number generators. +-- /Communications of the ACM/ 46(5):90–93. +-- +-- +-- * Thomas, D.B.; Leong, P.G.W.; Luk, W.; Villasenor, J.D. +-- (2007). Gaussian random number generators. +-- /ACM Computing Surveys/ 39(4). +-- diff --git a/Statistics/Resampling.hs b/Statistics/Resampling.hs --- a/Statistics/Resampling.hs +++ b/Statistics/Resampling.hs @@ -17,12 +17,12 @@ ) where import Control.Monad (forM_) -import Control.Monad.ST (unsafeSTToIO) +import Control.Monad.ST (ST) import Data.Array.Vector import Data.Array.Vector.Algorithms.Intro (sort) import Statistics.Function (createU) +import Statistics.RandomVariate (Gen, uniform) import Statistics.Types (Estimator, Sample) -import System.Random.Mersenne (MTGen, random) -- | A resample drawn randomly, with replacement, from a set of data -- points. Distinct from a normal array to make it harder for your @@ -33,20 +33,19 @@ -- | Resample a data set repeatedly, with replacement, computing each -- estimate over the resampled data. -resample :: MTGen -> [Estimator] -> Int -> Sample -> IO [Resample] +resample :: Gen s -> [Estimator] -> Int -> Sample -> ST s [Resample] resample gen ests numResamples samples = do - results <- unsafeSTToIO . mapM (const (newMU numResamples)) $ ests + results <- mapM (const (newMU numResamples)) $ ests loop 0 (zip ests results) - unsafeSTToIO $ do - mapM_ sort results - mapM (fmap Resample . unsafeFreezeAllMU) results + mapM_ sort results + mapM (fmap Resample . unsafeFreezeAllMU) results where loop k ers | k >= numResamples = return () | otherwise = do re <- createU n $ \_ -> do - r <- random gen + r <- uniform gen return (indexU samples (abs r `mod` n)) - unsafeSTToIO . forM_ ers $ \(est,arr) -> + forM_ ers $ \(est,arr) -> writeMU arr k . est $ re loop (k+1) ers n = lengthU samples diff --git a/Statistics/Resampling/Bootstrap.hs b/Statistics/Resampling/Bootstrap.hs --- a/Statistics/Resampling/Bootstrap.hs +++ b/Statistics/Resampling/Bootstrap.hs @@ -20,7 +20,7 @@ import Control.Exception (assert) import Data.Array.Vector (foldlU, filterU, indexU, lengthU) import Statistics.Distribution.Normal -import Statistics.Distribution (cumulative, inverse) +import Statistics.Distribution (cumulative, quantile) import Statistics.Resampling (Resample(..), jackknife) import Statistics.Sample (mean) import Statistics.Types (Estimator, Sample) @@ -75,9 +75,9 @@ hi = min (cumn a2) (ni - 1) where a2 = bias + b2 / (1 - accel * b2) b2 = bias - z1 - z1 = inverse standard ((1 - confidenceLevel) / 2) + z1 = quantile standard ((1 - confidenceLevel) / 2) cumn = round . (*n) . cumulative standard - bias = inverse standard (probN / n) + bias = quantile standard (probN / n) where probN = fromIntegral . lengthU . filterU ( Sample -> Double centralMoment a xs | a < 0 = error "Statistics.Sample.centralMoment: negative input" @@ -125,18 +129,21 @@ -- its mass is on the right of the distribution, with the tail on the -- left. -- --- > skewness . powers 3 $ toU [1,100,101,102,103] +-- > skewness $ toU [1,100,101,102,103] -- > ==> -1.497681449918257 -- -- A sample with positive skew is said to be /right-skewed/. -- --- > skewness . powers 5 $ toU [1,2,3,4,100] +-- > skewness $ toU [1,2,3,4,100] -- > ==> 1.4975367033335198 -- -- A sample's skewness is not defined if its 'variance' is zero. -- -- This function performs two passes over the sample, so is not subject -- to stream fusion. +-- +-- For samples containing many values very close to the mean, this +-- function is subject to inaccuracy due to catastrophic cancellation. skewness :: Sample -> Double skewness xs = c3 * c2 ** (-1.5) where c3 :*: c2 = centralMoments 3 2 xs @@ -152,6 +159,9 @@ -- -- This function performs two passes over the sample, so is not subject -- to stream fusion. +-- +-- For samples containing many values very close to the mean, this +-- function is subject to inaccuracy due to catastrophic cancellation. kurtosis :: Sample -> Double kurtosis xs = c4 / (c2 * c2) - 3 where c4 :*: c2 = centralMoments 4 2 xs diff --git a/Statistics/Sample/Powers.hs b/Statistics/Sample/Powers.hs --- a/Statistics/Sample/Powers.hs +++ b/Statistics/Sample/Powers.hs @@ -206,7 +206,8 @@ -- -- * Besset, D.H. (2000) Elements of statistics. -- /Object-oriented implementation of numerical methods/ --- pp. 311–331. +-- ch. 9, pp. 311–331. +-- -- -- * Anderson, G. (2009) Compute /k/th central moments in one -- pass. /quantblog/. diff --git a/statistics.cabal b/statistics.cabal --- a/statistics.cabal +++ b/statistics.cabal @@ -1,5 +1,5 @@ name: statistics -version: 0.2.2 +version: 0.3 synopsis: A library of statistical types, data, and functions description: This library provides a number of common functions and types useful @@ -15,6 +15,8 @@ . Computing with sample data: quantile estimation, kernel density estimation, bootstrap methods, and autocorrelation analysis. + . + Random variate generation under several different distributions. license: BSD3 license-file: LICENSE homepage: http://darcs.serpentine.com/statistics @@ -42,6 +44,7 @@ Statistics.KernelDensity Statistics.Math Statistics.Quantile + Statistics.RandomVariate Statistics.Resampling Statistics.Resampling.Bootstrap Statistics.Sample @@ -52,7 +55,7 @@ build-depends: base < 5, erf, - mersenne-random, + time, uvector >= 0.1.0.4, uvector-algorithms >= 0.2 if impl(ghc >= 6.10) @@ -62,6 +65,6 @@ -- gather extensive profiling data for now ghc-prof-options: -auto-all - ghc-options: -Wall -funbox-strict-fields -O2 + ghc-options: -Wall -funbox-strict-fields -O2 -fvia-C if impl(ghc >= 6.8) ghc-options: -fwarn-tabs