diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2014, Alp Mestanogullari + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Alp Mestanogullari nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,89 @@ +probable +======== + +Simple random value generation for haskell, using an efficient +random generator and minimizing system calls. But the library also +lets you work with distributions over a finite set, adapting +code from Eric Kidd's posts, and all the usual distributions +covered in the [statistics](http://hackage.haskell.org/package/statistics) +package. + +You can see how it looks in [examples](https://github.com/alpmestan/probable/tree/master/examples), or below. You can view the documentation for 0.1 [here](http://alpmestan.com/probable/). + +## Example + +Simple example of random generation for your types, using _probable_. + +``` haskell +module Main where + +import Control.Applicative +import Control.Monad +import Math.Probable + +import qualified Data.Vector.Unboxed as VU + +data Person = Person + { age :: Int + , weight :: Double + , salary :: Int + } deriving (Eq, Show) + +person :: RandT IO Person +person = + Person <$> intIn (1, 100) + <*> doubleIn (2, 130) + <*> intIn (500, 10000) + +randomPersons :: Int -> IO [Person] +randomPersons n = mwc $ listOf n person + +randomDoubles :: Int -> IO (VU.Vector Double) +randomDoubles n = mwc $ vectorOf n double + +main :: IO () +main = do + randomPersons 10 >>= mapM_ print + randomDoubles 10 >>= VU.mapM_ print +``` + +Distributions over finite sets, conditional probabilities and random sampling. + +``` haskell +module Main where + +import Math.Probable + +import qualified Data.Vector as V + +data Book = Interesting + | Boring + deriving (Eq, Show) + +bookPrior :: Finite d => d Book +bookPrior = weighted [ (Interesting, 0.2) + , (Boring, 0.8) + ] + +twoBooks :: Finite d => d (Book, Book) +twoBooks = do + book1 <- bookPrior + book2 <- bookPrior + return (book1, book2) + +sampleBooks :: RandT IO (V.Vector Book) +sampleBooks = vectorOf 10 bookPrior + +oneInteresting :: Fin (Book, Book) +oneInteresting = bayes $ do + (b1, b2) <- twoBooks + condition (b1 == Interesting || b2 == Interesting) + return (b1, b2) + +main :: IO () +main = do + print $ exact bookPrior + mwc sampleBooks >>= print + print $ exact twoBooks + print $ exact oneInteresting +``` diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/bench/distrib.hs b/bench/distrib.hs new file mode 100644 --- /dev/null +++ b/bench/distrib.hs @@ -0,0 +1,55 @@ +module Main where + +import Control.Applicative +import Control.Monad +import Control.Monad.Primitive +import Criterion.Main +import System.Random.MWC + +import qualified Data.Vector as V +import qualified Data.Vector.Unboxed as U +import Math.Probable + +n :: Int +n = 100000 + +data Fruit = Apple | Banana | Orange + deriving (Eq, Show) + +fruitDist :: FromFinite d => d Fruit +fruitDist = weighted [ (Apple, 0.3) + , (Banana, 0.6) + , (Orange, 0.1) + ] +liftF' :: Fin a -> IO a +liftF' dist = do + d <- withSystemRandom . asGenIO $ uniform + pick (P d) (exact dist) + +instance FromFinite IO where + weighted = liftF' . weighted + +observe1 :: RandT IO Fruit -> IO (V.Vector Fruit) +observe1 = mwc . vectorOf5 n + +observe2 :: RandT IO Fruit -> IO (V.Vector Fruit) +observe2 r = mwc . RandT $ \gen -> + V.replicateM n (runRandT r gen) + +observe4 :: IO Fruit + -> IO (V.Vector Fruit) +observe4 fr = + V.replicateM n (fr gen) + + +-- | Time to benchmark! +main :: IO () +main = do + defaultMain + [ + bgroup "big vector of fruits" + [ bench "observe1" $ whnfIO (observe1 fruitDist) + , bench "observe2" $ whnfIO (observe2 fruitDist) + , bench "observe3" $ whnfIO (observe4 fruitDist) + ] + ] diff --git a/bench/random.hs b/bench/random.hs new file mode 100644 --- /dev/null +++ b/bench/random.hs @@ -0,0 +1,88 @@ +module Main where + +import Control.Applicative +import Control.Monad +import Criterion.Main +import System.Random.MWC + +import qualified System.Random.MWC.Monad as MWCMonad + +import qualified Data.Vector as V +import qualified Data.Vector.Unboxed as U +import Math.Probable + +probable3 :: (U.Unbox a, Variate a) => Int -> IO (U.Vector a) +probable3 n = mwc (vectorOf3 n) +{-# INLINE probable3 #-} + +probable5 :: U.Unbox a => Int -> RandT IO a -> IO (U.Vector a) +probable5 n gen = mwc (vectorOf5 n gen) +{-# INLINE probable5 #-} + +mwc1 :: (U.Unbox a, Variate a) => Int -> IO (U.Vector a) +mwc1 n = + withSystemRandom . asGenIO $ + \gen -> uniformVector gen n +{-# INLINE mwc1 #-} + +mwc2 :: (U.Unbox a, Variate a) => Int -> IO (U.Vector a) +mwc2 n = + withSystemRandom . asGenIO $ + \gen -> U.replicateM n (uniform gen) +{-# INLINE mwc2 #-} + +mwcm :: (U.Unbox a, Variate a) => Int -> IO (U.Vector a) +mwcm n = + MWCMonad.runWithSystemRandom . MWCMonad.asRandIO $ + MWCMonad.toRand (flip uniformVector n) +{-# INLINE mwcm #-} + +{- +-- | Dummy 'Person' type +data Person = Person + { age :: Int + , weight :: Double + , salary :: Int + } deriving (Eq, Show) + +person :: (Generator g m Double, Generator g m Int) + => RandT g m Person +person = + Person <$> sampleUniform (1, 100) + <*> sampleUniform (2, 130) + <*> sampleUniform (500, 10000) + +randomPersons :: Int -> IO (V.Vector Person) +randomPersons n = mwc (vectorOf n person) +-} + +n :: Int +n = 10000000 + +-- | Time to benchmark! +main :: IO () +main = do + defaultMain + [ + bgroup "big vector of int" + [ bench "probable3" $ whnfIO (i $ probable3 n) + , bench "probable5" $ whnfIO (probable5 n int) + , bench "mwc-random" $ whnfIO (i $ mwc1 n) + , bench "mwc-random2" $ whnfIO (i $ mwc2 n) + , bench "mwc-random-monad2" $ whnfIO (i $ mwcm n) + ], + + bgroup "big vector of double" + [ bench "probable3" $ whnfIO (d $ probable3 n) + , bench "probable5" $ whnfIO (probable5 n double) + , bench "mwc-random" $ whnfIO (d $ mwc1 n) + , bench "mwc-random2" $ whnfIO (d $ mwc2 n) + , bench "mwc-random-monad2" $ whnfIO (d $ mwcm n) + ] + ] + +d :: IO (U.Vector Double) -> IO (U.Vector Double) +d = id + +i :: IO (U.Vector Int) -> IO (U.Vector Int) +i = id diff --git a/examples/finite.hs b/examples/finite.hs new file mode 100644 --- /dev/null +++ b/examples/finite.hs @@ -0,0 +1,36 @@ +module Main where + +import Math.Probable + +import qualified Data.Vector as V + +data Book = Interesting + | Boring + deriving (Eq, Show) + +bookPrior :: Finite d => d Book +bookPrior = weighted [ (Interesting, 0.2) + , (Boring, 0.8) + ] + +twoBooks :: Finite d => d (Book, Book) +twoBooks = do + book1 <- bookPrior + book2 <- bookPrior + return (book1, book2) + +sampleBooks :: RandT IO (V.Vector Book) +sampleBooks = vectorOf 10 bookPrior + +oneInteresting :: Fin (Book, Book) +oneInteresting = bayes $ do + (b1, b2) <- twoBooks + condition (b1 == Interesting || b2 == Interesting) + return (b1, b2) + +main :: IO () +main = do + print $ exact bookPrior + mwc sampleBooks >>= print + print $ exact twoBooks + print $ exact oneInteresting diff --git a/examples/montyhall.hs b/examples/montyhall.hs new file mode 100644 --- /dev/null +++ b/examples/montyhall.hs @@ -0,0 +1,127 @@ +{-# LANGUAGE BangPatterns #-} +module Main where + +import Data.List +import Math.Probable + +-- | We have 3 distinct doors +data Door = D1 | D2 | D3 + deriving (Eq, Show) + +-- | Our 3 doors in a list +doors :: [Door] +doors = [D1, D2, D3] + +-- | We can either Win or Lose +data Result = Win | Lose + deriving (Eq, Show) + +-- | Knowing the door we have chosen so far +-- and the one that's been opened +-- we decide to keep the one we've picked +keep :: Door -> Door -> Door +keep chosen _ = chosen + +-- | Knowing the door we have chosen so far +-- and the one that's been opened +-- we decide to switch to the third one +switch :: Door -> Door -> Door +switch chosen opened = head $ doors \\ [chosen, opened] + +-- | Given one of the two functions above ("strategies"), +-- and (the chosen door, the already opened one, and +-- the one with the car) +-- we check wether we won or not +resultOf :: (Door -> Door -> Door) -> Door -> Door -> Door -> Result +resultOf strategy chosen opened cardoor = + case strategy chosen opened == cardoor of + True -> Win + False -> Lose + +-- | We run the distribution of results and group the 'Win's and the 'Lose's +-- respective probabilities together +-- maybe this should be in the library, with a more general type... +collect :: Fin Result -> (Event Result, Event Result) +collect = f . exact + where f = toPair . foldl' combine (0, 0) + combine (!winP, !loseP) (Event Win p) = (winP+p, loseP) + combine (!winP, !loseP) (Event Lose p) = (winP, loseP+p) + toPair (winP, loseP) = (Event Win winP, Event Lose loseP) + +-- | Given a strategy to adopt, what's the distribution of Win/Lose ? +result :: (Door -> Door -> Door) + -> Fin Result +result strategy = do + -- we pick a door uniformly for hiding the car + carDoor <- uniformly doors + + -- we pick a door uniformly for the player + chosenDoor <- uniformly doors + + -- we open a door that neither hides the car (for suspense) + -- nor the one the player has picked (but they can be the same door) + openedDoor <- uniformly $ doors \\ [carDoor, chosenDoor] + + -- ok, the player tells us whether he decides to keep or switch + -- we now check whether he wins or not + let res = resultOf strategy chosenDoor openedDoor carDoor + + -- we return the result, to make this a distribution of results + return res + +-- | Given a strategy to adopt, distribution of doors for a Win ? +-- this uses Bayes' rule, and consequently lives in 'FinBayes' +result' :: (Door -> Door -> Door) + -> FinBayes (Door, Door, Door) +result' strategy = do + -- we pick a door uniformly for hiding the car + carDoor <- uniformly doors + + -- we pick a door uniformly for the player + chosenDoor <- uniformly doors + + -- we open a door that neither hides the car (for suspense) + -- nor the one the player has picked (but they can be the same door) + openedDoor <- uniformly $ doors \\ [carDoor, chosenDoor] + + -- ok, the player tells us whether he decides to keep or switch + -- we now check whether he wins or not + let res = resultOf strategy chosenDoor openedDoor carDoor + + -- here we discard all the combinations that don't lead to Win + condition (res == Win) + + -- and return the combination + return (chosenDoor, openedDoor, carDoor) + +main :: IO () +main = do + putStrLn $ "Using the conservative strategy: " + ++ show (collect $ result keep) + -- Using the conservative strategy: (Event Win 33.3%,Event Lose 66.7%) + + putStrLn $ "Switching: " + ++ show (collect $ result switch) + -- Switching: (Event Win 66.7%,Event Lose 33.3%) + + putStrLn "---" + + putStrLn $ "Winning (initial door, opened door, car door)'s - CONSERVATIVE:" + mapM_ print . exact . bayes $ result' keep + -- Event (D1,D2,D1) 16.7% + -- Event (D1,D3,D1) 16.7% + -- Event (D2,D1,D2) 16.7% + -- Event (D2,D3,D2) 16.7% + -- Event (D3,D1,D3) 16.7% + -- Event (D3,D2,D3) 16.7% + + putStrLn "---" + + putStrLn $ "Winning (initial door, opened door, car door)'s - SWITCHING:" + mapM_ print . exact . bayes $ result' switch + -- Event (D2,D3,D1) 16.7% + -- Event (D3,D2,D1) 16.7% + -- Event (D1,D3,D2) 16.7% + -- Event (D3,D1,D2) 16.7% + -- Event (D1,D2,D3) 16.7% + -- Event (D2,D1,D3) 16.7% diff --git a/examples/simple.hs b/examples/simple.hs new file mode 100644 --- /dev/null +++ b/examples/simple.hs @@ -0,0 +1,30 @@ +module Main where + +import Control.Applicative +import Control.Monad +import Math.Probable + +import qualified Data.Vector.Unboxed as VU + +data Person = Person + { age :: Int + , weight :: Double + , salary :: Int + } deriving (Eq, Show) + +person :: RandT IO Person +person = + Person <$> uniformIn (1, 100) + <*> uniformIn (2, 130) + <*> uniformIn (500, 10000) + +randomPersons :: Int -> IO [Person] +randomPersons n = mwc $ listOf n person + +randomDoubles :: Int -> IO (VU.Vector Double) +randomDoubles n = mwc $ vectorOf n double + +main :: IO () +main = do + randomPersons 10 >>= mapM_ print + randomDoubles 10 >>= VU.mapM_ print diff --git a/probable.cabal b/probable.cabal new file mode 100644 --- /dev/null +++ b/probable.cabal @@ -0,0 +1,121 @@ +name: probable +version: 0.1.0.0 +synopsis: Easy and reasonably efficient probabilistic programming and random generation +description: Easy and reasonably efficient probabilistic programming and random generation + . + This library gives a common language to speak about + probability distributions and + random generation, by wrapping both, when necessary, + in a 'RandT' monad defined in @Math.Probable.Random@. + This module also provides a lot of useful little + combinators for easily describing how random values for your + types should be generated. + . + In @Math.Probable.Distribution@, you'll find functions for generating + random values that follow any distribution supported by + . + . + In @Math.Probable.Distribution.Finite@, you'll find an adaptation + of Eric Kidd's work on probability monads (from + ). + . + You may want to check the examples bundled with this package, + viewable online at . + One of these examples is simple enough to be worth reproducing here. + . + > module Main where + > + > import Control.Applicative + > import Control.Monad + > import Math.Probable + > + > import qualified Data.Vector.Unboxed as VU + > + > data Person = Person Int -- ^ age + > Double -- ^ weight (kgs) + > Double -- ^ salary (e.g euros) + > deriving (Eq, Show) + > + > person :: RandT IO Person + > person = + > Person <$> uniformIn (1, 100) + > <*> uniformIn (2, 130) + > <*> uniformIn (500, 10000) + > + > randomPersons :: Int -> IO [Person] + > randomPersons n = mwc $ listOf n person + > + > randomDoubles :: Int -> IO (VU.Vector Double) + > randomDoubles n = mwc $ vectorOf n double + > + > main :: IO () + > main = do + > randomPersons 10 >>= mapM_ print + > randomDoubles 10 >>= VU.mapM_ print + . + Please report any feature request or problem, either by email + or through github's issues/feature requests. +homepage: http://github.com/alpmestan/probable +bug-reports: http://github.com/alpmestan/probable/issues +license: BSD3 +license-file: LICENSE +author: Alp Mestanogullari +maintainer: alpmestan@gmail.com +copyright: 2014 Alp Mestanogullari +category: Math, Statistics +build-type: Simple +cabal-version: >=1.10 +tested-with: GHC == 7.6.3, GHC == 7.8.2 +extra-source-files: bench/*.hs, + examples/*.hs, + README.md + + +source-repository head + type: git + location: https://github.com/alpmestan/probable.git + +library + exposed-modules: Math.Probable, + Math.Probable.Distribution, + Math.Probable.Distribution.Finite, + Math.Probable.Random + other-modules: + build-depends: base >=4.5 && <4.8, + statistics >= 0.10, + vector >= 0.10, + mwc-random >= 0.10, + primitive, + mtl, + transformers >= 0.3 + + hs-source-dirs: src + default-language: Haskell2010 + ghc-options: -O2 -funbox-strict-fields -Wall + +benchmark random + main-is: random.hs + hs-source-dirs: bench + ghc-options: -O2 -funbox-strict-fields + type: exitcode-stdio-1.0 + build-depends: base >= 4 && < 5, + vector >= 0.7, + criterion, + probable, + mwc-random, + mwc-random-monad + default-language: Haskell2010 + +benchmark distrib + main-is: distrib.hs + hs-source-dirs: bench + ghc-options: -O2 -funbox-strict-fields + type: exitcode-stdio-1.0 + build-depends: base >= 4 && < 5, + vector >= 0.7, + criterion, + probable, + mwc-random, + primitive + + default-language: Haskell2010 diff --git a/src/Math/Probable.hs b/src/Math/Probable.hs new file mode 100644 --- /dev/null +++ b/src/Math/Probable.hs @@ -0,0 +1,23 @@ +-- | +-- Module : Math.Probable +-- Copyright : (c) 2014 Alp Mestanogullari +-- License : BSD3 +-- Maintainer : alpmestan@gmail.com +-- Stability : experimental +-- Portability : GHC +-- +-- Easy, composable and efficient random number generation, +-- with support for distributions over a finite set +-- and your usual probability distributions. + +module Math.Probable + ( -- * random value generation + module Math.Probable.Random + + , -- * distributions + module Math.Probable.Distribution + ) where + +import Math.Probable.Distribution +import Math.Probable.Random + diff --git a/src/Math/Probable/Distribution.hs b/src/Math/Probable/Distribution.hs new file mode 100644 --- /dev/null +++ b/src/Math/Probable/Distribution.hs @@ -0,0 +1,282 @@ +module Math.Probable.Distribution + ( -- * Common distributions + beta + , cauchy + , cauchyStd + , chiSquared + , fisher + , gamma + , improperGamma + , geometric + , geometric0 + , student + , uniform + , normal + , standard + , normalFromSample + , exponential + , exponentialFromSample + + -- * Finite distributions + , module Math.Probable.Distribution.Finite + + -- * Utility functions + , continuous + , discrete + + ) where + +import Control.Monad.Primitive +import Math.Probable.Distribution.Finite +import Math.Probable.Random + +import Statistics.Distribution (ContGen, DiscreteGen, genContVar, genDiscreteVar) +import Statistics.Distribution.Beta +import Statistics.Distribution.CauchyLorentz +import qualified Statistics.Distribution.ChiSquared as Chi +import qualified Statistics.Distribution.Exponential as E +import Statistics.Distribution.FDistribution +import Statistics.Distribution.Gamma +import qualified Statistics.Distribution.Geometric as G +import qualified Statistics.Distribution.Normal as N +import Statistics.Distribution.StudentT +import Statistics.Distribution.Uniform +import Statistics.Types (Sample) + +-- | Sample from a continuous distribution from the 'statistics' package +-- +-- > λ> import qualified Statistics.Distribution.Normal as Normal +-- > λ> mwc $ continuous (Normal.normalDistr 0 1) +-- > -0.7266583064693862 +-- +-- This is equivalent to using 'normal' from this module. +continuous :: (ContGen d, PrimMonad m) + => d -- ^ the continuous distribution to sample from + -> RandT m Double +continuous d = RandT $ genContVar d +{-# INLINE continuous #-} + +-- | Sample from a discrete distribution from the 'statistics' package +-- +-- > λ> import qualified Statistics.Distribution.Normal as Normal +-- > λ> mwc $ discrete (Geo.geometric 0.6) +-- > 2 +-- +-- This is equivalent to using 'geometric' from this module. +discrete :: (DiscreteGen d, PrimMonad m) + => d -- ^ the discrete distribution to sample from + -> RandT m Int +discrete d = RandT $ genDiscreteVar d + +-- | Beta distribution (from @Statistics.Distribution.Beta@) +-- +-- > λ> mwc $ listOf 10 (beta 81 219) +-- > [ 0.23238372272745833,0.252972980515086,0.22708315774257903 +-- > , 0.25807200295967214,0.29794072226119983,0.24534701159196015 +-- > , 0.24766870269839578,0.2994199351220346,0.2728157476212405,0.2593318159573564 +-- > ] +beta :: PrimMonad m + => Double -- ^ shape parameter alpha + -> Double -- ^ shape parameter beta + -> RandT m Double +beta alpha bet = continuous $ betaDistr alpha bet + +-- | Cauchy distribution (from @Statistics.Distribution.Cauchy@) +-- +-- > λ> mwc $ listOf 10 (cauchy 0 0.1) +-- > [ -0.3932758718373347,0.490467375093784,4.2620417667423555e-2 +-- > , 3.370509874905657e-2,-8.186484692937862e-2,9.371858212168262e-2 +-- > , -1.1095818809115384e-2,3.0353983716155386e-2,0.22759697862410477 +-- > , -0.1881828277028582 ] +cauchy :: PrimMonad m + => Double -- ^ central point + -> Double -- ^ scale parameter + -> RandT m Double +cauchy p l = continuous $ cauchyDistribution p l + +-- | Cauchy distribution around 0, with scale 1 (from @Statistics.Distribution.Cauchy@) +-- +-- > λ> mwc $ listOf 10 cauchyStd +-- > [ 9.409701589649838,-7.361963972107541,0.168746305673769 +-- > , 5.091825420838711,-0.326080163135388,-1.2989850787629456 +-- > , -2.685658063444485,0.22671438734899435,-1.602349559644217e-2 +-- > , -0.6476292643908057 ] +cauchyStd :: PrimMonad m + => RandT m Double +cauchyStd = cauchy 0 1 + +-- | Chi-squared distribution (from @Statistics.Distribution.ChiSquared@) +-- +-- > λ> mwc $ listOf 10 (chiSquare 4) +-- > [ 8.068852054279787,1.861584389294606,6.3049415103095265 +-- > , 1.0512164068833838,1.6243237867165086,5.284901049954076 +-- > , 0.4773242487947021,1.1753876666374887,5.21554771873363 +-- > , 3.477574639460651 ] +chiSquared :: PrimMonad m + => Int -- ^ number of degrees of freedom + -> RandT m Double +chiSquared = continuous . Chi.chiSquared + +-- | Fisher's F-Distribution (from @Statistics.Distribution.FDistribution@) +-- +-- > λ> mwc $ listOf 10 (fisher 4 3) +-- > [ 3.437898578540642,0.844120450719367,1.9907851466347173 +-- > , 2.0089975147012784,1.3729208790549117,0.9380430357924707 +-- > , 2.642389323945247,1.0918121624055352,0.45650856735477335 +-- > , 2.5134537326659196 ] +fisher :: PrimMonad m + => Int + -> Int + -> RandT m Double +fisher a b = continuous $ fDistribution a b + +-- | Gamma distribution (from @Statistics.Distribution.Gamma@) +-- +-- > λ> mwc $ listOf 10 (gamma 3 0.1) +-- > [ 5.683745415884202e-2,0.20726188766138176,0.3150672538487696 +-- > , 0.4250825346490057,0.5586516230326105,0.46897413151474315 +-- > , 0.18374916962208182,9.93000480494153e-2,0.6057279704154832 +-- > , 0.11070190282993911 ] +gamma :: PrimMonad m + => Double -- ^ shape parameter k + -> Double -- ^ scale parameter theta + -> RandT m Double +gamma k theta = continuous $ gammaDistr k theta + +-- | Gamma distribution, without checking whether the parameter are valid +-- (from @Statistics.Distribution.Gamma@) +-- +-- > λ> mwc $ listOf 10 (improperGamma 3 0.1) +-- > [ 0.30431838005485,0.4044380297376584,2.8950141419406657e-2 +-- > , 0.468271612850147,0.18587792578128381,0.22735854572527045 +-- > , 0.5168050216325927,5.896911236207261e-2,0.24654560998405564 +-- > , 0.10557650513145429 ] +improperGamma :: PrimMonad m + => Double -- ^ shape parameter k + -> Double -- ^ scale parameter theta + -> RandT m Double +improperGamma k theta = continuous $ improperGammaDistr k theta + +-- | Geometric distribution. +-- +-- Distribution of the number of trials needed to get one success. +-- See @Statistics.Distribution.Geometric@ +-- +-- > λ> mwc $ listOf 10 (geometric 0.8) +-- > [2,1,1,1,1,1,1,2,1,5] +geometric :: PrimMonad m + => Double -- ^ success rate + -> RandT m Int +geometric = discrete . G.geometric + +-- | Geometric distribution. +-- +-- Distribution of the number of failures before getting one success. +-- See @Statistics.Distribution.Geometric@ +-- +-- > λ> mwc $ listOf 10 (geometric0 0.8) +-- > [0,0,0,0,0,1,1,0,0,0] +geometric0 :: PrimMonad m + => Double + -> RandT m Int +geometric0 = discrete . G.geometric0 + +-- | Student-T distribution (from @Statistics.Distribution.StudentT@) +-- +-- > λ> mwc $ listOf 10 (student 0.2) +-- > [ -14.221373473810829,-29.395749168822267,19.448665112984997 +-- > , -30.00446058929083,-0.5033202547957609,2.321975597874013 +-- > , 0.7884787761643617,-0.1895113832448149,-131.12901170537924 +-- > , 1.371956948317759 ] +student :: PrimMonad m + => Double + -> RandT m Double +student = continuous . studentT + +-- | Uniform distribution between 'a' and 'b' (from @Statistics.Distribution.Uniform@) +-- +-- > λ> mwc $ listOf 10 (uniform 0.1 0.2) +-- > [ 0.1711914559256124,0.1275212181343327,0.15347702635758945 +-- > , 0.1743662387063698,0.12047749686635312,0.10719840237585587 +-- > , 0.10543681342025846,0.13482973080648325,0.19779298960413577 +-- > , 0.1681037592576508 ] +uniform :: PrimMonad m + => Double + -> Double + -> RandT m Double +uniform a b = continuous $ uniformDistr a b + +-- | Normal distribution (from @Statistics.Distribution.Normal@) +-- +-- > λ> mwc $ listOf 10 (normal 4 1) +-- > [ 3.6815394812555144,3.5958531529526727,3.775960990625964 +-- > , 4.413109650155896,4.825826384709198,4.805629590118984 +-- > , 5.259267547365003,4.45410634165052,4.886537243027636 +-- > , 3.0409409067356954 ] +normal :: PrimMonad m + => Double -- ^ mean + -> Double -- ^ standard deviation + -> RandT m Double +normal mean stddev = continuous $ N.normalDistr mean stddev + +-- | The standard normal distribution (mean = 0, stddev = 1) (from @Statistics.Distribution.Normal@) +-- +-- > λ> mwc $ listOf 10 standard +-- > [ 0.2252627935262769,1.1831885443897947,-0.6577353418647461 +-- > , 2.1574536855051853,-0.16983072710637676,0.9667954287638821 +-- > , -1.8758605246293683,-0.8578048838241616,1.9516838769731923 +-- > , 0.43752574431460434 ] +standard :: PrimMonad m + => RandT m Double +standard = continuous N.standard + +-- | Create a normal distribution using parameters estimated from the sample +-- (from @Statistics.Distribution.Normal@) +-- +-- > λ> mwc . listOf 10 $ +-- > normalFromSample $ +-- > V.fromList [1,1,1,3,3,3,4 +-- > ,4,4,4,4,4,4,4 +-- > ,4,4,4,4,4,5,5 +-- > ,5,7,7,7] +-- > [ 7.1837511677441395,2.388433817342809,5.252282321156134 +-- > , 4.988163140851522,0.40102386713467864,4.4840751065620665 +-- > , 2.1471370686776874,2.6591948802201046,3.843667372514598 +-- > , 1.7650436484843248 ] +normalFromSample :: PrimMonad m + => Sample -- ^ sample + -> RandT m Double +normalFromSample = continuous . N.normalFromSample + +-- | Exponential distribution (from @Statistics.Distribution.Exponential@) +-- +-- > λ> mwc $ listOf 10 (exponential 0.2) +-- > [ 5.713524665694821,1.7774315204594584,2.434017573227628 +-- > , 5.463202731505528,0.5403008025455847,14.346316301765576 +-- > , 7.380393612391503,24.800854500680032,0.8731076703020924 +-- > , 6.1661076502236645 ] +exponential :: PrimMonad m + => Double -- ^ lambda (scale) parameter + -> RandT m Double +exponential = continuous . E.exponential + +-- | Exponential distribution given a sample (from @Statistics.Distribution.Exponential@) +-- +-- > λ> mwc $ listOf 10 (exponentialFromSample $ V.fromList [1,1,1,0]) +-- > [ 0.4237050903604833,1.934301502525168,0.7435728843566659 +-- > , 1.8720263209574293,0.605750265970631,0.24103955067365979 +-- > , 0.6294952762436511,1.660404952631443,0.6448230847113577 +-- > , 0.8891555734786789 ] +exponentialFromSample :: PrimMonad m + => Sample + -> RandT m Double +exponentialFromSample = continuous . E.exponentialFromSample + + + + + + + + + diff --git a/src/Math/Probable/Distribution/Finite.hs b/src/Math/Probable/Distribution/Finite.hs new file mode 100644 --- /dev/null +++ b/src/Math/Probable/Distribution/Finite.hs @@ -0,0 +1,385 @@ +{-# LANGUAGE GeneralizedNewtypeDeriving, + BangPatterns, + TupleSections, + FlexibleInstances, + TypeSynonymInstances #-} +-- | +-- Module : Math.Probable +-- License : BSD3 +-- Maintainer : alpmestan@gmail.com +-- Stability : experimental +-- Portability : GHC +-- +-- /Fun with finite distributions!/ +-- +-- This all pretty much comes from Eric Kidd's series +-- of blog posts at . +-- +-- I have adapted it a bit by making it fit into my own +-- random generation/sampling scheme. +-- +-- The idea and purpose of this module should be clear after going +-- through an example. First, let's import the library and 'vector'. +-- +-- > import Math.Probable +-- > import qualified Data.Vector as V +-- +-- We are going to talk about Books, and particularly about whether a given +-- book is interesting or not. +-- +-- > data Book = Interesting +-- > | Boring +-- > deriving (Eq, Show) +-- +-- Let's say we have very particular tastes, and that we think +-- that only 20% of all books are interesting (that's not so small actually. oh well). +-- +-- > bookPrior :: Finite d => d Book +-- > bookPrior = weighted [ (Interesting, 0.2) +-- > , (Boring, 0.8) +-- > ] +-- +-- 'weighted' belongs to the 'Finite' class, which represents +-- types that can somehow represent a distribution over a finite set. +-- That makes our distribution polymorphic in how we will use it. Awesome! +-- +-- So how does it look? +-- +-- > λ> exact bookPrior -- in ghci +-- > [Event Interesting 20.0%,Event Boring 80.0%] +-- +-- 'exact' takes 'Fin' 'a' and gives you the +-- inner list that 'Fin' uses to represent the distribution. +-- +-- Now, what if we pick two books? First, how do we even do that? +-- Well, any instance of 'Finite' must be a 'Monad', so you have your +-- good old /do notation/. The ones provided by this package also +-- provide 'Functor' and 'Applicative' instances, but let's use +-- do. +-- +-- > twoBooks :: Finite d => d (Book, Book) +-- > twoBooks = do +-- > book1 <- bookPrior +-- > book2 <- bookPrior +-- > return (book1, book2) +-- +-- Nothing impressive. We pick a book with the prior +-- we defined above, then another, pair them together +-- and hand the pair back. What this will actually do +-- is behave just like in the list monad, but in addition +-- to this it will combine the probabilities of the various +-- events we could be dealing with in the appropriate way. +-- +-- So, how about we verify what I just said: +-- +-- > λ> exact twoBooks +-- > [ Event (Interesting,Interesting) 4.0% +-- > , Event (Interesting,Boring) 16.0% +-- > , Event (Boring,Interesting) 16.0% +-- > , Event (Boring,Boring) 64.0% +-- > ] +-- +-- Nice! Let's take a look at a more complicated scenario now. +-- +-- What if we wanted to take a look at the same distribution, +-- with just a difference: we want at least one of the books to +-- be an Interesting one. +-- +-- > oneInteresting :: Fin (Book, Book) +-- > oneInteresting = bayes $ do -- notice the call to bayes +-- > (b1, b2) <- twoBooks +-- > condition (b1 == Interesting || b2 == Interesting) +-- > return (b1, b2) +-- +-- We get two books from the previous distribution, and use 'condition' +-- to restrict the current distribution to the values of b1 and b2 +-- that verify our condition. This lifts us in the 'FinBayes' type, +-- where our probabilistic computations can "fail" in some sense. +-- If you want to discard values and restrict the ones on which you'll +-- run further computations, use 'condition'. +-- +-- However, how do we view the distribution now, without having all +-- those 'Maybe's in the middle? That's what 'bayes' is for. It runs +-- the computations for the distribution and discards all the ones +-- where any 'condition' wasn't satisfied. In particular, it means +-- it hands you back a normal 'Fin' distribution. +-- +-- If we run this one: +-- +-- > λ> exact oneInteresting +-- > [ Event (Interesting,Interesting) 11.1% +-- > , Event (Interesting,Boring) 44.4% +-- > , Event (Boring,Interesting) 44.4% +-- > ] +-- +-- Note that these finite distribution types support random sampling too: +-- +-- * If one of your distributions has a type like "Finite d => d X", +-- you can actually consider it as a 'RandT' value, from which you can sample. +-- +-- * If you have a 'Fin' distribution, you can use 'liftF' (lift 'Fin') +-- to randomly sample an element from it, by more or less following +-- the distribution's probabilities. +-- +-- > -- example of the former +-- > sampleBooks :: RandT IO (V.Vector Book) +-- > sampleBooks = vectorOf 10 bookPrior +-- +-- > λ> mwc sampleBooks +-- > fromList [Interesting,Boring,Boring,Boring,Boring +-- > ,Boring,Boring,Interesting,Boring,Boring] +-- +-- > λ> mwc $ listOf 4 (liftF oneInteresting) -- example of the latter +-- > [ (Boring,Interesting) +-- > , (Boring,Interesting) +-- > , (Boring,Interesting) +-- > , (Interesting,Boring) +-- > ] +module Math.Probable.Distribution.Finite + ( -- * Probability type + P(..), prob + + -- * 'Event' type + , Event(..), never + -- * 'EventT' monad transformer + , EventT(..) + -- * Finite distributions: 'Finite' and 'Fin' + , Finite(..), Fin, exact, uniformly, liftF + -- * Bayes' rule: 'FinBayes' + , FinBayes, bayes, condition, onlyJust + ) where + +import Control.Applicative +import Control.Monad +import Control.Monad.Primitive +import Control.Monad.Trans +import Control.Monad.Trans.Maybe +import Data.List +import Math.Probable.Random + +-- | Probability type: wrapper around Double +-- for a nicer Show instance and for more easily +-- enforcing normalization of weights +newtype P = P Double + deriving (Eq, Ord, Fractional, Num, Real, RealFrac) + +-- | Get the underlying probability +-- +-- > λ> prob (P 0.1) +-- > 0.1 +prob :: P -> Double +prob (P x) = x + +instance Show P where + show (P p) = show intPart ++ "." ++ show fracPart ++ "%" + where digits = round (1000 * p) + intPart = digits `div` 10 :: Int + fracPart = digits `mod` 10 :: Int + +-- | An event, and its probability +data Event a = Event a {-# UNPACK #-} !P + deriving (Eq, Show) + +-- | This event never happens (probability of 0) +-- +-- > never = Event undefined 0 +never :: Event a +never = Event undefined 0 + +instance Functor Event where + fmap f (Event evt p) = Event (f evt) p + +instance Applicative Event where + pure evt = Event evt 1 + + Event f p1 <*> Event e p2 + | p1 == 0 || p2 == 0 = never + | otherwise = Event (f e) (p1*p2) + +instance Monad Event where + return evt = Event evt 1 + + (Event evt p) >>= f | p == 0 = never + | otherwise = Event e' (p*p') + where Event e' p' = f evt + + +-- | 'EventT' monad transformer +-- +-- It pairs a value with a probability within the 'm' monad +newtype EventT m a = EventT { runEventT :: m (Event a) } + +instance Monad m => Functor (EventT m) where + fmap = liftM + +instance (Functor m, Monad m) => Applicative (EventT m) where + pure = return + + mf <*> x = mf >>= \f -> fmap f x + +instance Monad m => Monad (EventT m) where + return = lift . return + + m >>= f = EventT go + where go = do + ph <- runEventT m + case ph of + Event e p1 | p1 == 0 -> return never + | otherwise -> + do Event e' p2 <- runEventT (f e) + return $ Event e' (p1 * p2) + + +instance MonadTrans EventT where + lift x = EventT (liftM return x) + +-- | Create a 'Finite' distribution over the values in +-- the list, each with an equal probability +-- +-- > λ> exact $ uniformly [True, False] +-- > [Event True 50.0%,Event False 50.0%] +uniformly :: Finite d => [a] -> d a +uniformly = weighted . map (,1) +{-# INLINE uniformly #-} + +-- | 'Fin' is just 'EventT []' +-- +-- You can think of 'Fin a' meaning '[Event a]' +-- i.e a list of the possible outcomes of type 'a' +-- with their respective probability +type Fin = EventT [] + +-- | See the outcomes of a finite distribution and their probabilities +-- +-- > λ> exact $ uniformly [True, False] +-- > [Event True 50.0%,Event False 50.0%] +-- +-- > λ> data Fruit = Apple | Orange deriving (Eq, Show) +-- > λ> exact $ uniformly [Apple, Orange] +-- > [Event Apple 50.0%,Event Orange 50.0%] +-- +-- > λ> exact $ weighted [(Apple, 0.8), (Orange, 0.2)] +-- > [Event Apple 80.0%,Event Orange 20.0%] +exact :: Fin a -> [Event a] +exact = runEventT +{-# INLINE exact #-} + +-- | 'FinBayes' is 'Fin' with a 'MaybeT' layer +-- +-- What is that for? The 'MaybeT' lets us express +-- the fact that what we've drawn from the distribution +-- isn't of interest anymore, using 'condition', +-- and observing the remaining cases, using 'bayes', +-- to get back to a normal finite distribution. Example: +-- +-- > data Wine = Good | Bad deriving (Eq, Show) +-- > +-- > wines :: Finite d => d Wine +-- > wines = weighted [(Good, 0.2), (Bad, 0.8)] +-- > +-- > twoWines :: Finite d => d (Wine, Wine) +-- > twoWines = (,) <*> wines <$> wines +-- > +-- > decentMeal :: FinBayes (Wine, Wine) +-- > decentMeal = do +-- > (wine1, wine2) <- twoWines +-- > -- we only consider the outcomes of 'twoWines' +-- > -- where at least one of the two wines is good +-- > -- because we're having a nice meal and are looking +-- > -- for a decent pair of wine +-- > condition (wine1 == Good || wine2 == Good) +-- > return (wine1, wine2) +-- > +-- > -- to view the distribution, applying +-- > -- Bayes' rule on our way: +-- > exact (bayes decentMeal) +type FinBayes = MaybeT Fin + +-- | This is the core of 'FinBayes'. If the 'Bool' is false, +-- the current computation is shortcuited (sent to a 'Nothing' +-- in 'MaybeT') and won't be included when running the distribution +-- with 'bayes'. See the documentation of 'FinBayes' for an example. +condition :: Bool -> FinBayes () +condition = MaybeT . return . toMaybe + where toMaybe True = Just () + toMaybe False = Nothing + +-- | This functions discards all the elements of the distribution +-- for which the call to 'condition' yielded 'Nothing'. +-- While 'condition' does the mapping to 'Maybe' values, +-- this function discards all of those values for which the condition +-- was not met. +bayes :: FinBayes a -> Fin a +bayes = onlyJust . runMaybeT + +-- | Keeps only the 'Just's and remove the 'Maybe' layer +-- in the distribution. +onlyJust :: Fin (Maybe a) -> Fin a +onlyJust dist + | total > 0 = EventT (map adjust filtered) + | otherwise = EventT [] + where filtered = catMaybes' (runEventT dist) + total = sum (map proba filtered) + adjust (Event x p) = Event x (p / total) + proba (Event _ p) = p + -- value (Event x _) = x + +-- | This function, used by 'onlyJust', discards all the events +-- holding a 'Nothing'. +catMaybes' :: [Event (Maybe a)] -> [Event a] +catMaybes' [] = [] +catMaybes' (Event Nothing _ : xs) = + catMaybes' xs +catMaybes' (Event (Just x) p : xs) = + Event x p : catMaybes' xs + +-- | T distribution of probabilities +-- over a finite set. +class (Functor d, Monad d) => Finite d where + + -- | The only requirement is to somehow + -- be able to represent the distribution + -- corresponding to the list given as argument, e.g: + -- + -- > weighted [(True, 0.8), (False, 0.2)] + -- + -- It should also be able to handle the normalization for you. + -- + -- > weighted [(True, 8), (False, 2)] + weighted :: [(a, Double)] -> d a + +instance Finite Fin where + weighted l = EventT $ map weight l + + where weight (x, w) = Event x $ P (w / total) + total = foldl' (\w (_, w') -> w + w') 0 l + +instance PrimMonad m => Finite (RandT m) where + weighted = liftF . weighted + +-- | Make finite distributions ('Fin') citizens of +-- 'RandT' by simply sampling an element at random +-- while still approximately preserving the distribution +-- +-- > λ> mwc . liftF $ uniformly [True, False] +-- > False +-- > λ> mwc . liftF $ uniformly [True, False] +-- > True +-- > λ> mwc . liftF $ weighted [("Haskell", 99), ("PHP", 1)] +-- > "Haskell" +liftF :: PrimMonad m => Fin a -> RandT m a +liftF dist = do + d <- double + pick (P d) (exact dist) + +-- | Look for an 'Event' in that list that has +-- a probability superior to the one given as an +-- argument. +pick :: Monad m => P -> [Event a] -> m a +pick _ [] = error "pick: no value to pick" +pick n (Event x p : evts) + | n <= p = return x + | otherwise = pick (n-p) evts + +instance Finite FinBayes where + weighted = lift . weighted diff --git a/src/Math/Probable/Random.hs b/src/Math/Probable/Random.hs new file mode 100644 --- /dev/null +++ b/src/Math/Probable/Random.hs @@ -0,0 +1,483 @@ +{-# LANGUAGE BangPatterns, + TypeFamilies #-} + +-- | +-- Module : Math.Probable +-- Copyright : (c) 2014 Alp Mestanogullari +-- License : BSD3 +-- Maintainer : alpmestan@gmail.com +-- Stability : experimental +-- Portability : GHC +-- +-- Random number generation based on 'MWC.Gen', +-- defined as a Monad transformer. +-- +-- Quickstart, in ghci: +-- +-- > λ> import Math.Probable +-- > λ> import Control.Applicative +-- > λ> mwc double +-- > 0.2756820707828763 +-- > λ> mwc word64 +-- > 12175918187293541909 +-- > λ> mwc $ (,) <$> bool <*> intIn (0, 10) +-- > (True,7) +-- > λ> mwc $ do { n <- intIn (1, 10) ; listOf n (listOf 2 bool) } +-- > [ [False,True],[True,False],[False,True],[False,False],[False,False], +-- > [False,True],[True,False],[True,False],[True,True],[False,False] +-- > ] +-- +-- This module features a bunch of combinators that can help you create +-- some random generation descriptions easily, and in a very familiar style. +-- +-- You can easily combine them through the 'Monad' instance for 'RandT' +-- which really just make sure everyone gets a 'MWC.Gen' (from mwc-random) +-- eventually. This of course makes 'RandT' a 'Functor' and an 'Applicative'. +-- +-- > import Math.Probable +-- > +-- > data Person = +-- > Person { name :: String +-- > , age :: Int +-- > , salary :: Double +-- > } +-- > deriving (Eq, Show) +-- > +-- > randomPerson :: PrimMonad m +-- > => RandT m Person +-- > randomPerson = do +-- > -- we pick a random length +-- > -- for the person's name +-- > nameLen <- intIn (3, 10) +-- > +-- > -- and just express what a random Person +-- > -- should be, Applicative-style +-- > Person <$> pickName nameLen -- pick a name +-- > <*> intIn (0, 100) -- an Int between 0 and 100 +-- > <*> doubleIn (0, 10000) -- a Double between 0 and 10000 +-- > +-- > where pickName nameLen = do +-- > -- the initial, between 'A' and 'Z' +-- > initial <- chr `fmap` intIn (65, 90) +-- > +-- > (initial:) `fmap` +-- > -- the rest, between 'a' and 'z' +-- > listOf (nameLen - 1) +-- > (chr `fmap` intIn (97, 122)) +-- +-- This is all nice, but how do we actually sample such a Person? +-- You just have to call 'mwc': +-- +-- > λ> mwc randomPerson +-- > Person {name = "Ojeesra", age = 83, salary = 3075.9945184521885} +-- +-- So any value of type 'RandT m a' is something that you'll eventually +-- run in 'm' (hence 'IO' or 'ST' 's') for generating a /random value/ of +-- type 'a'. Note that 'mwc' forces the execution using 'withSystemRandom' +-- and gets you back in 'IO', whereas 'mwcST' gets you back in 'ST' 's'. +-- +-- My simple name generation routine can help you pick a name for your baby, +-- if you are having one soon. +-- +-- > λ> map name `fmap` mwc (listOf 10 randomPerson) +-- > ["Npujbc","Faidx","Zusha","Ghbipic","Ljaestei","Fktcfonnxe","Hlvkolds","Zpws","Zgmrkrdv","Rhcd"] +-- +-- If we were to make a generator that could generate more familiar +-- and creativity-free names, we wouldn't sample uniformly +-- from the alphabet. + +module Math.Probable.Random + ( -- * 'RandT' type + RandT(..) + + , -- * Actually generating random values + mwc, mwcST + + , -- * Combinators for generating individual values + uniformIn + + , int, int8, int16, int32, int64 + , intIn, int8In, int16In, int32In, int64In + + , word, word8, word16, word32, word64 + , wordIn, word8In, word16In, word32In, word64In + + , float, double + , floatIn, doubleIn + + , bool + + , -- * Filling containers with random values + listOf, vectorOf, vectorOfVariate + ) where + +import Control.Applicative +import Control.Monad.Identity +import Control.Monad.Primitive +import Control.Monad.ST +import Data.Int +import Data.Word +import qualified Data.Vector.Generic as G +import qualified System.Random.MWC as MWC + +-- | 'RandT' type, equivalent to a +-- 'ReaderT' ('MWC.Gen' ('PrimState' m)) +-- +-- This lets you build simple or complex random generation +-- routines without having the generator passed all around +-- and just run the whole thing in the end, most likely +-- by using 'mwc'. +newtype RandT m a = + RandT { runRandT :: MWC.Gen (PrimState m) -> m a } + +instance Monad m => Monad (RandT m) where + return x = RandT $ \_ -> return x + {-# INLINE return #-} + + (RandT g) >>= f = + RandT $ \gen -> do + !v <- g gen + !res <- runRandT (f v) gen + return res + {-# INLINE (>>=) #-} + +instance Monad m => Functor (RandT m) where + fmap f r = RandT $ \gen -> return . f =<< runRandT r gen + {-# INLINE fmap #-} + +instance Monad m => Applicative (RandT m) where + pure = return + {-# INLINE pure #-} + + (<*>) = ap + {-# INLINE (<*>) #-} + +-- | Generate a random 'Int'. The whole 'Int' range is used. +-- +-- > λ> mwc int +-- > 8354496680947360541 +int :: PrimMonad m => RandT m Int +int = RandT MWC.uniform +{-# INLINE int #-} + +-- | Generate a random 'Int' in the given range. +-- +-- > λ> mwc $ intIn (0, 10) +-- > 7 +intIn :: PrimMonad m => (Int, Int) -> RandT m Int +intIn (a, b) = RandT $ MWC.uniformR (a, b) +{-# INLINE intIn #-} + +-- | Generate a random 'Int8'. The whole 'Int8' range is used. +-- +-- > λ> mwc int8 +-- > -65 +int8 :: PrimMonad m => RandT m Int8 +int8 = RandT MWC.uniform +{-# INLINE int8 #-} + +-- | Generate a random 'Int8' in the given range +-- +-- > λ> mwc $ int8In (-10, 10) +-- > -3 +int8In :: PrimMonad m => (Int8, Int8) -> RandT m Int8 +int8In (a, b) = RandT $ MWC.uniformR (a, b) +{-# INLINE int8In #-} + +-- | Generate a random 'Int16'. The whole 'Int16' range is used. +-- +-- > λ> mwc int16 +-- > 15413 +int16 :: PrimMonad m => RandT m Int16 +int16 = RandT MWC.uniform +{-# INLINE int16 #-} + +-- | Generate a random 'Int16' in the given range +-- +-- > λ> mwc $ int16In (-500, 30129) +-- > 9501 +int16In :: PrimMonad m => (Int16, Int16) -> RandT m Int16 +int16In (a, b) = RandT $ MWC.uniformR (a, b) +{-# INLINE int16In #-} + +-- | Generate a random 'Int32'. The whole 'Int32' range is used. +-- +-- > λ> mwc int32 +-- > 1774441747 +int32 :: PrimMonad m => RandT m Int32 +int32 = RandT MWC.uniform +{-# INLINE int32 #-} + +-- | Generate a random 'Int32' in the given range. +-- +-- > λ> mwc $ int32In (-500, 30129) +-- > 8012 +int32In :: PrimMonad m => (Int32, Int32) -> RandT m Int32 +int32In (a, b) = RandT $ MWC.uniformR (a, b) +{-# INLINE int32In #-} + +-- | Generate a random 'Int64'. The whole 'Int64' range is used. +-- +-- > λ> mwc int64 +-- > -2596387699802756017 +int64 :: PrimMonad m => RandT m Int64 +int64 = RandT MWC.uniform +{-# INLINE int64 #-} + +-- | Generate a random 'Int64' in the given range. +-- +-- > λ> mwc $ int64In (-2^30, 30) +-- > -630614786 +int64In :: PrimMonad m => (Int64, Int64) -> RandT m Int64 +int64In (a, b) = RandT $ MWC.uniformR (a, b) +{-# INLINE int64In #-} + +-- | Generate a random 'Word'. The whole 'Word' range is used. +-- +-- > λ> mwc word +-- > 3106215968599504888 +word :: PrimMonad m => RandT m Word +word = RandT MWC.uniform +{-# INLINE word #-} + +-- | Generate a random 'Word' in the given range. +-- +-- > λ> mwc $ wordIn (1, 64) +-- > 28 +wordIn :: PrimMonad m => (Word, Word) -> RandT m Word +wordIn (a, b) = RandT $ MWC.uniformR (a, b) +{-# INLINE wordIn #-} + +-- | Generate a random 'Word8'. The whole 'Word8' range is used. +-- +-- > λ> mwc word8 +-- > 231 +word8 :: PrimMonad m => RandT m Word8 +word8 = RandT MWC.uniform +{-# INLINE word8 #-} + +-- | Generate a random 'Word8' in the given range +-- +-- > λ> mwc $ word8In (2, 15) +-- > 3 +word8In :: PrimMonad m => (Word8, Word8) -> RandT m Word8 +word8In (a, b) = RandT $ MWC.uniformR (a, b) +{-# INLINE word8In #-} + +-- | Generate a random 'Word16'. The whole 'Word16' range is used. +-- +-- > λ> mwc word16 +-- > 31127 +word16 :: PrimMonad m => RandT m Word16 +word16 = RandT MWC.uniform +{-# INLINE word16 #-} + +-- | Generate a random 'Word16' in the given range. +-- +-- > λ> mwc $ word16In (2^13, 2^14) +-- > 8885 +word16In :: PrimMonad m => (Word16, Word16) -> RandT m Word16 +word16In (a, b) = RandT $ MWC.uniformR (a, b) +{-# INLINE word16In #-} + +-- | Generate a random 'Word32'. The whole 'Word32' range is used. +-- +-- > λ> mwc word32 +-- > 3917666696 +word32 :: PrimMonad m => RandT m Word32 +word32 = RandT MWC.uniform +{-# INLINE word32 #-} + +-- | Generate a random 'Word32' in the given range. +-- +-- > λ> mwc $ word32In (100, 330) +-- > 125 +word32In :: PrimMonad m => (Word32, Word32) -> RandT m Word32 +word32In (a, b) = RandT $ MWC.uniformR (a, b) +{-# INLINE word32In #-} + +-- | Generate a random 'Word64'. The whole 'Word64' range is used. +-- +-- > λ> mwc word64 +-- > 12496697905424132339 +word64 :: PrimMonad m => RandT m Word64 +word64 = RandT MWC.uniform +{-# INLINE word64 #-} + +-- | Generate a random 'Word64' in the given range. +-- +-- > λ> mwc $ word64In (2^45, 2^46) +-- > 59226619151303 +word64In :: PrimMonad m => (Word64, Word64) -> RandT m Word64 +word64In (a, b) = RandT $ MWC.uniformR (a, b) +{-# INLINE word64In #-} + +-- | Generate a random 'Float' between 0 (excluded) +-- and 1 (included) +-- +-- > λ> mwc float +-- > 0.11831179 +float :: PrimMonad m => RandT m Float +float = RandT MWC.uniform +{-# INLINE float #-} + +-- | Generate a random 'Float' in the given range +-- +-- > λ> mwc $ floatIn (0.20, 3.14) +-- > 1.3784513 +floatIn :: PrimMonad m => (Float, Float) -> RandT m Float +floatIn (a, b) = RandT $ MWC.uniformR (a, b) +{-# INLINE floatIn #-} + +-- | Generate a random 'Double' between 0 (excluded) +-- and 1 (included) +-- +-- > λ> mwc double +-- > 0.7689412928620208 +double :: PrimMonad m => RandT m Double +double = RandT MWC.uniform +{-# INLINE double #-} + +-- | Generate a random 'Double' in the given range +-- +-- > λ> mwc $ doubleIn (-30.121121445, 0.129898878612) +-- > -13.612464813256999 +doubleIn :: PrimMonad m => (Double, Double) -> RandT m Double +doubleIn (a, b) = RandT $ MWC.uniformR (a, b) +{-# INLINE doubleIn #-} + +-- | Generate a random 'Bool' +-- +-- > λ> mwc bool +-- > False +bool :: PrimMonad m => RandT m Bool +bool = RandT MWC.uniform +{-# INLINE bool #-} + +-- | A generic function for sampling uniformly any type +-- that implements 'MWC.Variate'. +-- +-- All the 'xxxIn' functions from this module just +-- call 'MWC.uniformR'. +uniformIn :: (MWC.Variate a, PrimMonad m) => (a, a) -> RandT m a +uniformIn (a, b) = RandT $ MWC.uniformR (a, b) +{-# INLINE uniformIn #-} + +-- | Take a 'RandT' value and run it in 'IO', +-- generating all the random values described by +-- the 'RandT'. It just uses 'MWC.withSystemRandom' +-- so you really should try hard to put your whole +-- random generation logic in 'RandT' and call +-- 'mwc' in the end, thus initialising the generator +-- only once and generating everything with it. +-- +-- See the documentation for 'MWC.withSystemRandom' for more about this. +-- +-- > λ> mwc $ (+2) `fmap` int8 +-- > 34 +mwc :: RandT IO a + -> IO a +mwc = MWC.withSystemRandom + . MWC.asGenIO + . runRandT +{-# INLINE mwc #-} + +-- | If for some reason you have a 'RandT' ('ST' 's') +-- you can run it from 'IO' just like we do with 'mwc'. +-- +-- > λ> mwcST $ listOf 4 bool +-- > [False,False,True,True] +mwcST :: RandT (ST s) a + -> IO a +mwcST = MWC.withSystemRandom + . MWC.asGenST + . runRandT +{-# INLINE mwcST #-} + +-- | Repeatedly run a random computation +-- yielding a value of type 'a' to get +-- a list of random values of type 'a'. +-- +-- > λ> mwc (listOf 30 float) +-- > [ 5.438623e-2,0.78114086,0.4954672,0.5958733,0.47243807,5.883485e-2 +-- > , 5.500287e-2,0.79262286,0.5528683,0.7628807,0.80705905,0.15368962 +-- > , 0.8654971,0.4560417,0.23922172,0.5069659,0.8130155,0.6559351 +-- > , 1.31405e-2,0.25705606,0.7134138,0.79111993,0.7529769,0.10573909 +-- > , 0.37731406,0.6289338,0.85156864,0.15691182,0.9910314,8.133593e-2 +-- > ] +-- +-- > λ> mwc (sum `fmap` listOf 30 float) +-- > 15.037931 +listOf :: Monad m + => Int + -> RandT m a + -> RandT m [a] +listOf n r = replicateM n r +{-# INLINE listOf #-} + +-- | A function for generating a vector +-- of the given length for values +-- whose types are instances of 'MWC.Variate'. +-- +-- This function is generic in the type of vector it returns, +-- any instance of 'G.Vector' will do. +-- +-- It's just a wrapper arround 'MWC.uniformVector' +-- and doesn't really use the 'Monad' instance of 'RandT'. +-- +-- But if you want to have a vector of 'Person's, +-- you have to use 'vectorOf'. +-- +-- > λ> import qualified Data.Vector.Unboxed as V +-- > λ> :set -XScopedTypeVariables +-- > λ> v :: V.Vector Double <- mwc $ vectorOfVariate 10 +-- > λ> V.mapM_ print v +-- > 3.8565084196117705e-2 +-- > 0.575103826646098 +-- > 0.379710162825715 +-- > 0.4066991135077237 +-- > 0.9778431248247549 +-- > 0.3786223745680838 +-- > 0.4361789615081698 +-- > 0.9904407826187301 +-- > 0.2951087330670904 +-- > 0.1533350329892028 +vectorOfVariate :: (PrimMonad m, MWC.Variate a, G.Vector v a) + => Int + -> RandT m (v a) +vectorOfVariate n = + RandT $ \gen -> MWC.uniformVector gen n +{-# INLINE vectorOfVariate #-} + +-- | A function for generating a vector of the given +-- length with random values /of any type/ +-- (in contrast to 'vectorOfVariate'). +-- +-- It is generic in the 'G.Vector' instance it +-- hands you back. It's implemented in terms of +-- 'G.replicateM' and has been benchmarked to perform +-- as well as 'MWC.uniformVector' on simple types +-- ('MWC.uniformVector' can't generate values for types +-- that don't have a 'MWC.Variate' instance). +-- +-- > λ> import qualified Data.Vector.Unboxed as V +-- > λ> :set -XScopedTypeVariables +-- > λ> v :: V.Vector Int <- mwc $ vectorOf 10 int +-- > λ> V.mapM_ print v +-- > -3920053790769159788 +-- > 3983393642052845448 +-- > 1528310798822685910 +-- > 3522283620461337684 +-- > 6451017362937898910 +-- > 1929485210691770214 +-- > 8547527164583329795 +-- > 3298785082692387491 +-- > 4019024417224980311 +-- > -5216301990322376953 +vectorOf :: (Monad m, G.Vector v a) + => Int + -> RandT m a + -> RandT m (v a) +vectorOf n r = + RandT $ \gen -> G.replicateM n (runRandT r gen) +{-# INLINE vectorOf #-}