diff --git a/System/Random/Effect.hs b/System/Random/Effect.hs deleted file mode 100644 --- a/System/Random/Effect.hs +++ /dev/null @@ -1,90 +0,0 @@ -{-# LANGUAGE DeriveDataTypeable #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE FlexibleContexts #-} --- | A random number effect, using a pure mersenne twister under --- the hood. This should be plug-and-play with any application --- making use of extensible effects. --- --- Patches, even for the smallest of documentation bugs, are --- always welcome! -module System.Random.Effect ( Random - , mkRandom - , mkRandomIO - , randomInt - , randomInt64 - , randomWord - , randomWord64 - , randomDouble - ) where - -import Data.Int -import Data.Typeable -import Data.Word -import qualified System.Random.Mersenne.Pure64 as SR - -import Control.Eff -import Control.Eff.State -import Control.Eff.Lift - --- | A pure mersenne twister pseudo-random number generator. -data Random = Random {-# UNPACK #-} !SR.PureMT - deriving Typeable - --- | Create a random number generator from a 'Word64' seed. -mkRandom :: Word64 -> Random -mkRandom = Random . SR.pureMT -{-# INLINE mkRandom #-} - --- | Create a new random number generator, using the clocktime as the base for --- the seed. This must be called from a computation with a lifted base effect --- of 'IO'. -mkRandomIO :: Member (Lift IO) r - => Eff r Random -mkRandomIO = lift (fmap Random SR.newPureMT) -{-# INLINE mkRandomIO #-} - --- | Runs an effectful random computation, returning the computation's result. -runRandomState :: Random - -> Eff (State Random :> r) w - -> Eff r w -runRandomState seed computation = - fmap snd (runState seed computation) -{-# INLINE runRandomState #-} - --- | A generalized form of generating a random number of the correct type --- from System.Random.Mersenne.Pure64. -randomF :: Member (State Random) r - => (SR.PureMT -> (a, SR.PureMT)) - -> Eff r a -randomF f = do - (Random old) <- getState - let (val, new) = f old - putState (Random new) - return val -{-# INLINE randomF #-} - --- | Yield a new 'Int' value from the generator. The full 64 bits will be used --- on a 64 bit machine. -randomInt :: Member (State Random) r => Eff r Int -randomInt = randomF SR.randomInt -{-# INLINE randomInt #-} - --- | Yield a new 'Word' value from the generator. -randomWord :: Member (State Random) r => Eff r Word -randomWord = randomF SR.randomWord -{-# INLINE randomWord #-} - --- | Yield a new 'Int64' value from the generator. -randomInt64 :: Member (State Random) r => Eff r Int64 -randomInt64 = randomF SR.randomInt64 -{-# INLINE randomInt64 #-} - --- | Yield a new 'Word64' value from the generator. -randomWord64 :: Member (State Random) r => Eff r Word64 -randomWord64 = randomF SR.randomWord64 -{-# INLINE randomWord64 #-} - --- | Yield a new 53-bit precise 'Double' value from the generator. -randomDouble :: Member (State Random) r => Eff r Double -randomDouble = randomF SR.randomDouble -{-# INLINE randomDouble #-} diff --git a/src/System/Random/Effect.hs b/src/System/Random/Effect.hs new file mode 100644 --- /dev/null +++ b/src/System/Random/Effect.hs @@ -0,0 +1,540 @@ +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} +-- | A random number effect, using a pure mersenne twister under +-- the hood. This should be plug-and-play with any application +-- making use of extensible effects. +-- +-- Patches, even for the smallest of documentation bugs, are +-- always welcome! +module System.Random.Effect ( Random + -- * Seeding + , mkRandom + , mkRandomIO + -- * Running + , runRandomState + -- * Uniform Distributions + , uniformIntDist + , uniformRealDist + -- * Bernoulli Distributions + , bernoulliDist + , binomialDist + , negativeBinomialDist + , geometricDist + -- * Poisson Distributions + , poissonDist + , exponentialDist + , gammaDist + , weibullDist + , extremeValueDist + -- * Normal Distributions + , normalDist + , lognormalDist + , chiSquaredDist + , cauchyDist + , fisherFDist + , studentTDist + -- * Sampling Distributions + , DiscreteDistributionHelper + , buildDDH + , discreteDist + , piecewiseConstantDist + -- * Raw Generators + , randomInt + , randomInt64 + , randomWord + , randomWord64 + , randomDouble + , randomBits + , randomBitList + ) where + +import Control.Applicative +import Control.Monad +import Control.Monad.ST +import Data.Bits +import Data.Int +import Data.List +import Data.Ratio +import Data.Typeable +import Data.Vector.Algorithms.Search +import Data.Vector ( Vector ) +import qualified Data.Vector as V +import Data.Word +import Statistics.Distribution + +import qualified Statistics.Distribution.CauchyLorentz as DCL +import qualified Statistics.Distribution.ChiSquared as DC +import qualified Statistics.Distribution.Exponential as DE +import qualified Statistics.Distribution.FDistribution as DF +import qualified Statistics.Distribution.Gamma as DG +import qualified Statistics.Distribution.Normal as DN +--import qualified Statistics.Distribution.Poisson as DP +import qualified Statistics.Distribution.StudentT as DS + +import qualified System.Random.Mersenne.Pure64 as SR + +import Control.Eff +import Control.Eff.State +import Control.Eff.Lift + +-- | A pure mersenne twister pseudo-random number generator. +data Random = Random {-# UNPACK #-} !SR.PureMT + deriving Typeable + +-- | Create a random number generator from a 'Word64' seed. +mkRandom :: Word64 -> Random +mkRandom = Random . SR.pureMT +{-# INLINE mkRandom #-} + +-- | Create a new random number generator, using the clocktime as the base for +-- the seed. This must be called from a computation with a lifted base effect +-- of 'IO'. +mkRandomIO :: Member (Lift IO) r + => Eff r Random +mkRandomIO = lift (fmap Random SR.newPureMT) +{-# INLINE mkRandomIO #-} + +-- | Runs an effectful random computation, returning the computation's result. +runRandomState :: Random + -> Eff (State Random :> r) w + -> Eff r w +runRandomState seed computation = + fmap snd (runState seed computation) +{-# INLINE runRandomState #-} + +-- | A generalized form of generating a random number of the correct type +-- from System.Random.Mersenne.Pure64. +randomF :: Member (State Random) r + => (SR.PureMT -> (a, SR.PureMT)) + -> Eff r a +randomF f = do + (Random old) <- getState + let (val, new) = f old + putState (Random new) + return val +{-# INLINE randomF #-} + +-- | Yield a new 'Int' value from the generator. The full 64 bits will be used +-- on a 64 bit machine. +randomInt :: Member (State Random) r => Eff r Int +randomInt = randomF SR.randomInt +{-# INLINE randomInt #-} + +-- | Yield a new 'Word' value from the generator. +randomWord :: Member (State Random) r => Eff r Word +randomWord = randomF SR.randomWord +{-# INLINE randomWord #-} + +-- | Yield a new 'Int64' value from the generator. +randomInt64 :: Member (State Random) r => Eff r Int64 +randomInt64 = randomF SR.randomInt64 +{-# INLINE randomInt64 #-} + +-- | Yield a new 'Word64' value from the generator. +randomWord64 :: Member (State Random) r => Eff r Word64 +randomWord64 = randomF SR.randomWord64 +{-# INLINE randomWord64 #-} + +-- | Yield a new 53-bit precise 'Double' value from the generator. +-- The returned number will be in the range [0, 1). +randomDouble :: Member (State Random) r => Eff r Double +randomDouble = randomF SR.randomDouble +{-# INLINE randomDouble #-} + +-- | Yields a set of random from the internal generator, +-- using 'randomWord64' internally. +randomBits :: (Bits x, Member (State Random) r) + => Eff r x +randomBits = do + let z = clearBit (bit 0) 0 -- zero, so we can get the number of bits + nBits = bitSize z + + -- we OR with zero to get the number of bits above. + -- it shouldn't affect the output. + (z .|.) . bitsToInteger <$> randomBitList nBits +{-# INLINE randomBits #-} + +-- | Returns a list of bits which have been randomly generated. +randomBitList :: Member (State Random) r + => Int -- ^ The number of bits to generate + -> Eff r [Bool] +randomBitList k = do + let iters = (k `div` 64) + 1 + + breakBits w = map (testBit w) [0..(bitSize w - 1)] + + word64s <- replicateM iters randomWord64 + + return $ take k (concatMap breakBits word64s) + +-- | Returns the maximum set bit in an integer. +maxBit :: Integer -> Int +maxBit x = + let loop y !k + | y == 0 = k + | otherwise = loop (y `unsafeShiftR` 1) (k+1) + in loop x 0 + +-- | Repeat a computation until it succeeds a test. +loopUntil :: Monad m => (a -> Bool) -> m a -> m a +loopUntil f c = do + x <- c + case f x of + True -> return x + False -> loopUntil f c + +b2i :: Bits a => a -> Bool -> a +b2i x True = setBit x 0 +b2i x False = x +{-# INLINE b2i #-} + +bitsToInteger :: Bits a => [Bool] -> a +bitsToInteger = + let z = clearBit (bit 0) 0 + in foldl' (\x -> b2i (unsafeShiftL x 1)) z + +-- | Implementation of 'uniformIntDist', with shared calculations +-- with other uniformIntDists passed as parameters. This lets us +-- share work with other calls to 'uniformIntDist' with the same +-- parameters. +-- +-- Returns a number in the inclusive range [0, range]. +-- +-- 'hMask' is the mask to get the number of bits per random number +-- we generate. +-- +-- 'nBits' is the number of bits in 'range'. +uniformIntDist' :: Member (State Random) r + => Integer -- ^ range + -> Int -- ^ numBits + -> Eff r Integer +uniformIntDist' range nBits + | range == 0 = return 0 + | otherwise = + loopUntil (<= range) $ + bitsToInteger <$> randomBitList nBits + +-- | Generates a uniformly distributed random number in +-- the inclusive range [a, b]. +uniformIntDist :: Member (State Random) r + => Integer -- ^ a + -> Integer -- ^ b + -> Eff r Integer +uniformIntDist a' b' = + let a = min a' b' + b = max a' b' + range = b - a + maxB = maxBit range + in (a+) <$> uniformIntDist' range maxB +{-# INLINE uniformIntDist #-} + +-- | The part of 'uniformRealDist' that does all the work. +-- We factor it out so we can inline 'uniformRealDist', +-- and possibly share as much work as possible. +uniformRealDist' :: Member (State Random) r + => Double -- ^ a + -> Double -- ^ range + -> Eff r Double +uniformRealDist' a range = do + d <- randomDouble + return (d * range + a) + +-- | Generates a uniformly distributed random number in +-- the inclusive range [a, b]. +-- +-- NOTE: This code might not be correct, in that the +-- returned value may not be perfectly uniformly +-- distributed. If you know how to make one of +-- these a better way, PLEASE send me a pull request. +-- I just stole this implementation from the C++11 +-- header. +uniformRealDist :: Member (State Random) r + => Double -- ^ a + -> Double -- ^ b + -> Eff r Double +uniformRealDist a' b' = + let a = min a' b' + b = max a' b' + range = b - a + + in uniformRealDist' a range +{-# INLINE uniformRealDist #-} + +-- | Samples a continuous distribution. +-- +-- Generates random numbers as if they were sampled by the given +-- distribution. +-- +-- This is implemented with the inverse transform rule: +-- . +sampleContDist :: (Member (State Random) r + , ContDistr d) + => d -- ^ The distribution to sample. + -> Eff r Double +sampleContDist d = + quantile d <$> randomDouble + +-- | Produces random boolean values, according to a discrete probability. +-- +-- k must be in the range [0, 1]. +bernoulliDist :: Member (State Random) r + => Rational -- ^ k: The fraction of results which should be true. + -> Eff r Bool +bernoulliDist k + | k < 0 || k > 1 = + error $ "system-random-effect: bernoulliDist: fraction 'k' is out of range: " ++ show k + | k == 0 = return False + | otherwise = do + let (n, d) = (numerator k, denominator k) + dist <- uniformIntDist 1 d + return $ dist <= n +{-# INLINE bernoulliDist #-} + +-- | The value obtained is the number of successes in a sequence of +-- t yes/no experiments, each of which succeeds with probability p. +-- +-- t must be >= 0 +-- p must be in the range [0, 1]. +-- +-- Warning: NOT IMPLEMENTED! +binomialDist :: Member (State Random) r + => Int -- ^ t + -> Rational -- ^ p + -> Eff r Int +binomialDist t p + | p < 0 || p > 1 = + error $ "system-random-effect: binomialDist: fraction 'p' is out of range: " ++ show p + | t == 0 || p == 0 = return 0 + | p == 1 = return t + | otherwise = do + trials <- V.replicateM t randomDouble + + let p' = realToFrac p + succeeded = V.map (<= p') trials + numSuccesses = + V.foldl' (\s b -> if b then s+1 else s) 0 succeeded + + return numSuccesses + +-- | The value represents the number of failures in a series of +-- independent yes/no trials (each succeeds with probability p), +-- before exactly k successes occur. +-- +-- p must be in the range (0, 1] +-- k must be >= 0 +-- +-- Warning: NOT IMPLEMENTED! +negativeBinomialDist :: Member (State Random) r + => Rational -- ^ p + -> Integer -- ^ k + -> Eff r Integer +negativeBinomialDist p k + | p <= 0 || p > 1 = + error $ "system-random-effect: negativeBinomialDist: fraction 'p' is out of range: " ++ show p + | p == 1 || k == 0 = return 0 + | otherwise = + error "system-random-effect: negativeBinomialDist: TODO. Patches welcome!" + +-- | The value represents the number of yes/no trials (each +-- succeeding with probability p) which are necessary to +-- obtain a single success. +-- +-- 'geometricDist' p is equivalent to negativeBinomialDist 1 p +-- +-- p must be in the range (0, 1] +-- +-- Warning: NOT IMPLEMENTED! +geometricDist :: Member (State Random) r + => Rational -- ^ p + -> Eff r Integer +geometricDist p + | p <= 0 || p > 1 = + error $ "system-random-effect: geometricDist: fraction 'p' is out of range: " ++ show p + | p == 1 = return 0 + | otherwise = + error "system-random-effect: geometricDist: TODO: Patches welcome!" + +-- | The value obtained is the probability of exactly i +-- occurrences of a random event if the expected, mean +-- number of its occurrence under the same conditions +-- (on the same time/space interval) is μ. +-- +-- Warning: NOT IMPLEMENTED! +poissonDist :: Member (State Random) r + => Double -- ^ μ + -> Eff r Double -- ^ i +poissonDist = + error "system-random-effect: poissonDist: TODO: Patches welcome!" + +-- | The value obtained is the time/distance until the next +-- random event if random events occur at constant rate λ +-- per unit of time/distance. For example, this distribution +-- describes the time between the clicks of a Geiger counter +-- or the distance between point mutations in a DNA strand. +-- +-- This is the continuous counterpart of 'geometricDist'. +exponentialDist :: Member (State Random) r + => Double -- ^ λ. Scale parameter. + -> Eff r Double +exponentialDist lambda = + sampleContDist (DE.exponential lambda) + +-- | For floating-point α, the value obtained is the sum of α +-- independent exponentially distributed random variables, +-- each of which has a mean of β. +gammaDist :: Member (State Random) r + => Double -- ^ α. The shape parameter. + -> Double -- ^ β. The scale parameter. + -> Eff r Double +gammaDist alpha beta = + sampleContDist (DG.gammaDistr alpha beta) + +-- | ??? +-- +-- Warning: NOT IMPLEMENTED! +weibullDist :: Member (State Random) r + => Double -- ^ α. The shape parameter. + -> Double -- ^ β. The scale parameter. + -> Eff r Double +weibullDist = + error "system-random-effect: weibullDist: TODO: Patches welcome!" + +-- | ??? +-- +-- Warning: NOT IMPLEMENTED! +extremeValueDist :: Member (State Random) r + => Double -- ^ α. The shape parameter. + -> Double -- ^ β. The scale parameter. + -> Eff r Double +extremeValueDist = + error "system-random-effect: extremeValueDist: TODO: Patches welcome!" + +-- | Generates random numbers as sampled from the +-- normal distribution. +normalDist :: Member (State Random) r + => Double -- ^ μ. The mean. + -> Double -- ^ σ. The standard deviation. + -> Eff r Double +normalDist mu sigma = + sampleContDist (DN.normalDistr mu sigma) + +-- | Generates a log-normally distributed random number. +-- This is based off of sampling the normal distribution, +-- and then following the instructions at +-- . +lognormalDist :: Member (State Random) r + => Double -- ^ μ. The mean. + -> Double -- ^ σ. The standard deviation. + -> Eff r Double +lognormalDist mu sigma = do + z <- sampleContDist simpleND + return $ exp (mu + sigma*z) + +-- | Share this computation over every call of lognormalDist +simpleND :: DN.NormalDistribution +simpleND = DN.normalDistr 0 1 +{-# NOINLINE simpleND #-} + +-- | Produces random numbers according to a chi-squared distribution. +chiSquaredDist :: Member (State Random) r + => Int -- ^ n. The number of degrees of freedom. + -> Eff r Double +chiSquaredDist n + | n <= 0 = + error $ "system-random-effect: chiSquaredDist: invalid degrees of freedom: " ++ show n + | otherwise = + sampleContDist (DC.chiSquared n) + +-- | Produced random numbers according to a Cauchy (or Lorentz) distribution. +cauchyDist :: Member (State Random) r + => Double -- ^ Central point + -> Double -- ^ Scale parameter (full width half maximum) + -> Eff r Double +cauchyDist a b = + sampleContDist (DCL.cauchyDistribution a b) + +-- | Produces random numbers according to an F-distribution. +-- +-- m and n are the degrees of freedom. +fisherFDist :: Member (State Random) r + => Int -- ^ m + -> Int -- ^ n + -> Eff r Double +fisherFDist m n = + sampleContDist (DF.fDistribution m n) + +-- | This distribution is used when estimating the mean of an +-- unknown normally distributed value given n+1 independent +-- measurements, each with additive errors of unknown standard +-- deviation, as in physical measurements. Or, alternatively, +-- when estimating the unknown mean of a normal distribution +-- with unknown standard deviation, given n+1 samples. +studentTDist :: Member (State Random) r + => Double -- ^ The number of degrees of freedom + -> Eff r Double +studentTDist d = + sampleContDist (DS.studentT d) + +-- | Contains a sorted list of cumulative probabilities, so we +-- can do a sample by generating a uniformly distributed random +-- number in the range [0, 1), and binary searching the vector +-- for where to put it. +newtype DiscreteDistributionHelper = + DDH (Vector Rational) + +-- | Performs O(n) work building a table which we can later use +-- sample with 'discreteDist'. +buildDDH :: [Integer] -> DiscreteDistributionHelper +buildDDH xs = + let vs = V.fromList xs + s = V.sum vs + ns = V.map (% s) vs -- normalize the list + in DDH (V.postscanl' (+) 0 ns) + +-- | Given a pre-build 'DiscreteDistributionHelper' (use 'buildDDH'), +-- produces random integers on the interval [0, n), where the +-- probability of each individual integer i is defined as w_i/S, +-- that is the weight of the ith integer divided by the sum of all +-- n weights. +-- +-- i.e. This function produces an integer with probability equal to +-- the weight given in its index into the parameter to 'buildDDH'. +discreteDist :: Member (State Random) r + => DiscreteDistributionHelper + -> Eff r Int +discreteDist (DDH xs) = do + y <- realToFrac <$> randomDouble + return $ runST $ do + mv <- V.thaw xs + binarySearch mv y + +-- | This function produces random floating-point numbers, which +-- are uniformly distributed within each of the several subintervals +-- [b_i, b_(i+1)), each with its own weight w_i. The set of interval +-- boundaries and the set of weights are the parameters of this +-- distribution. +-- +-- For example, `piecewiseConstantDist [ 0, 1, 10, 15 ] +-- (buildDDH [ 1, 0, 1 ])` +-- will produce values between 0 and 1 half the time, and values +-- between 10 and 15 the other half of the time. +piecewiseConstantDist :: Member (State Random) r + => [Double] -- ^ Intervals + -> DiscreteDistributionHelper -- ^ Weights + -> Eff r Double +piecewiseConstantDist intervals weights@(DDH rs) + | V.length rs + 1 /= length intervals = + error $ "system-random-effect: piecewiseConstantDist:" + ++ " Incongruent parameter lengths." + ++ " intervals=" ++ show intervals + ++ " weights=" ++ show rs + | otherwise = do + idx <- discreteDist weights + + let vints = V.fromList intervals + (l, r) = (vints V.! idx, vints V.! (idx+1)) + + uniformRealDist l r diff --git a/system-random-effect.cabal b/system-random-effect.cabal --- a/system-random-effect.cabal +++ b/system-random-effect.cabal @@ -1,5 +1,5 @@ name: system-random-effect -version: 0.1.0.0 +version: 0.1.1.0 synopsis: Random number generation for extensible effects. homepage: https://github.com/wowus/system-random-effect license: BSD3 @@ -11,10 +11,37 @@ cabal-version: >=1.10 library + hs-source-dirs: src/ + ghc-options: -Wall exposed-modules: System.Random.Effect - build-depends: base >=4.6 && <4.7 + build-depends: base == 4.6.* , extensible-effects == 1.1.* , mersenne-random-pure64 == 0.2.* + , statistics == 0.10.* + , vector == 0.10.* + , vector-algorithms == 0.5.* default-language: Haskell2010 + +test-suite extensible-effects-tests + type: exitcode-stdio-1.0 + main-is: Test.hs + hs-source-dirs: test/ + + ghc-options: -rtsopts=all -threaded + + build-depends: base == 4.6.* + , QuickCheck == 2.* + , HUnit == 1.2.* + , test-framework == 0.8.* + , test-framework-hunit == 0.3.* + , test-framework-quickcheck2 == 0.3.* + , extensible-effects == 1.1.* + , system-random-effect + + default-language: Haskell2010 + +source-repository head + type: git + location: https://github.com/wowus/system-random-effect diff --git a/test/Test.hs b/test/Test.hs new file mode 100644 --- /dev/null +++ b/test/Test.hs @@ -0,0 +1,37 @@ +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE FlexibleContexts #-} +module Main ( main ) where + +import Control.Eff + +import System.Random.Effect + +import Control.Monad (void) +import Data.Word +import Data.Typeable + +import Test.Framework (defaultMain, testGroup) +import Test.Framework.Providers.HUnit +import Test.Framework.Providers.QuickCheck2 + +import Test.HUnit hiding (State) +import Test.QuickCheck + +main :: IO () +main = defaultMain tests + +checkRange :: (Integer, Integer) -> Integer -> Bool +checkRange (low, high) x = + x >= low && x <= high + +testUniformRandom :: Integer -> Integer -> Word64 -> Bool +testUniformRandom a b seed = + let low = min a b + high = max a b + + in checkRange (low, high) $ run $ runRandomState (mkRandom seed) $ do + uniformIntDist a b + +tests = + [ testProperty "random range" testUniformRandom + ]