mwc-probability

Sampling function-based probability distributions.
Log | Files | Refs | README | LICENSE

commit 6006519f4ded4a93de8bd9666a850d6d2fc52aab
Author: Jared Tobin <jared@jtobin.ca>
Date:   Sun, 25 May 2014 15:39:01 +1000

Initial commit.

Diffstat:
A.gitignore | 4++++
ALICENSE | 20++++++++++++++++++++
AREADME.md | 38++++++++++++++++++++++++++++++++++++++
ASetup.hs | 2++
Amwc-probability.cabal | 57+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/System/Random/MWC/Probability.hs | 116+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6 files changed, 237 insertions(+), 0 deletions(-)

diff --git a/.gitignore b/.gitignore @@ -0,0 +1,4 @@ +.cabal-sandbox +dist +*swp +cabal.sandbox.config diff --git a/LICENSE b/LICENSE @@ -0,0 +1,20 @@ +Copyright (c) 2014 Jared Tobin + +Permission is hereby granted, free of charge, to any person obtaining +a copy of this software and associated documentation files (the +"Software"), to deal in the Software without restriction, including +without limitation the rights to use, copy, modify, merge, publish, +distribute, sublicense, and/or sell copies of the Software, and to +permit persons to whom the Software is furnished to do so, subject to +the following conditions: + +The above copyright notice and this permission notice shall be included +in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY +CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, +TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE +SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. diff --git a/README.md b/README.md @@ -0,0 +1,38 @@ +# mwc-probability + +Sampling function-based probability distributions. + +A simple probability distribution type, where distributions are characterized +by sampling functions. + +This implementation is a thin layer over @mwc-random@, which handles RNG +state-passing automatically by using a @PrimMonad@ like @IO@ or @ST s@ under +the hood. + +Includes Functor, Applicative, Monad, and MonadTrans instances. + +Examples +-------- + +Transform a distribution's support while leaving its density structure +invariant: + + -- uniform over [0, 1] to uniform over [1, 2] + succ <$> uniform + +Sequence distributions together using bind: + + > -- a beta-binomial conjugate distribution + > beta 1 10 >>= binomial 10 + +Use do-notation to build complex joint distributions from composable, +local conditionals: + + > hierarchicalModel = do + > [c, d, e, f] <- replicateM 4 $ uniformR (1, 10) + > a <- gamma c d + > b <- gamma e f + > p <- beta a b + > n <- uniformR (5, 10) + > binomial n p + diff --git a/Setup.hs b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/mwc-probability.cabal b/mwc-probability.cabal @@ -0,0 +1,57 @@ +name: mwc-probability +version: 0.1.0.0 +synopsis: Sampling function-based probability distributions. +description: + + A simple probability distribution type, where distributions are characterized + by sampling functions. + . + This implementation is a thin layer over @mwc-random@, which handles RNG + state-passing automatically by using a @PrimMonad@ like @IO@ or @ST s@ under + the hood. + . + Includes Functor, Applicative, Monad, and MonadTrans instances. + . + /Examples/ + . + Transform a distribution's support while leaving its density structure + invariant: + . + > -- uniform over [0, 1] to uniform over [1, 2] + > succ <$> uniform + . + Sequence distributions together using bind: + . + > -- a beta-binomial conjugate distribution + > beta 1 10 >>= binomial 10 + . + Use do-notation to build complex joint distributions from composable, + local conditionals: + . + > hierarchicalModel = do + > [c, d, e, f] <- replicateM 4 $ uniformR (1, 10) + > a <- gamma c d + > b <- gamma e f + > p <- beta a b + > n <- uniformR (5, 10) + > binomial n p + +homepage: http://github.com/jtobin/mwc-probability +license: MIT +license-file: LICENSE +author: Jared Tobin +maintainer: jared@jtobin.ca +category: Math +build-type: Simple +cabal-version: >=1.10 + +library + exposed-modules: System.Random.MWC.Probability + default-language: Haskell2010 + hs-source-dirs: src + build-depends: + base >= 4.7 && < 4.8 + , mwc-random >= 0.13.1.1 + , mtl >= 2.2.0.1 + , primitive >= 0.5.3.0 + diff --git a/src/System/Random/MWC/Probability.hs b/src/System/Random/MWC/Probability.hs @@ -0,0 +1,116 @@ +{-# OPTIONS_GHC -Wall #-} + +module System.Random.MWC.Probability ( + module MWC + , Prob(..) + , uniform + , uniformR + , categorical + , standard + , normal + , logNormal + , exponential + , gamma + , inverseGamma + , chiSquare + , beta + , dirichlet + , bernoulli + , binomial + ) where + +import Control.Applicative +import Control.Monad +import Control.Monad.Primitive +import Control.Monad.Trans +import System.Random.MWC as MWC hiding (uniform, uniformR) +import qualified System.Random.MWC as QMWC +import qualified System.Random.MWC.Distributions as MWC.Dist + +newtype Prob m a = Prob { sample :: Gen (PrimState m) -> m a } + +instance Monad m => Functor (Prob m) where + fmap h (Prob f) = Prob $ liftM h . f + {-# INLINE fmap #-} + +instance Monad m => Applicative (Prob m) where + pure = return + (<*>) = ap + {-# INLINE pure #-} + {-# INLINE (<*>) #-} + +instance Monad m => Monad (Prob m) where + return = Prob . const . return + m >>= h = Prob $ \g -> do + z <- sample m g + sample (h z) g + {-# INLINE return #-} + {-# INLINE (>>=) #-} + +instance MonadTrans Prob where + lift m = Prob $ const m + {-# INLINE lift #-} + +uniform :: (PrimMonad m, Variate a) => Prob m a +uniform = Prob QMWC.uniform +{-# INLINE uniform #-} + +uniformR :: (PrimMonad m, Variate a) => (a, a) -> Prob m a +uniformR r = Prob $ QMWC.uniformR r +{-# INLINE uniformR #-} + +categorical :: PrimMonad m => [a] -> Prob m a +categorical cs = do + j <- uniformR (0, length cs - 1) + return $ cs !! j +{-# INLINE categorical #-} + +standard :: PrimMonad m => Prob m Double +standard = Prob MWC.Dist.standard +{-# INLINE standard #-} + +normal :: PrimMonad m => Double -> Double -> Prob m Double +normal m sd = Prob $ MWC.Dist.normal m sd +{-# INLINE normal #-} + +logNormal :: PrimMonad m => Double -> Double -> Prob m Double +logNormal m sd = exp <$> normal m sd +{-# INLINE logNormal #-} + +exponential :: PrimMonad m => Double -> Prob m Double +exponential r = Prob $ MWC.Dist.exponential r +{-# INLINE exponential #-} + +gamma :: PrimMonad m => Double -> Double -> Prob m Double +gamma a b = Prob $ MWC.Dist.gamma a b +{-# INLINE gamma #-} + +inverseGamma :: PrimMonad m => Double -> Double -> Prob m Double +inverseGamma a b = recip <$> gamma a b +{-# INLINE inverseGamma #-} + +chiSquare :: PrimMonad m => Int -> Prob m Double +chiSquare k = Prob $ MWC.Dist.chiSquare k +{-# INLINE chiSquare #-} + +beta :: PrimMonad m => Double -> Double -> Prob m Double +beta a b = do + u <- gamma a 1 + w <- gamma b 1 + return $ u / (u + w) +{-# INLINE beta #-} + +dirichlet :: PrimMonad m => [Double] -> Prob m [Double] +dirichlet as = do + zs <- mapM (`gamma` 1) as + return $ map (/ sum zs) zs +{-# INLINE dirichlet #-} + +bernoulli :: PrimMonad m => Double -> Prob m Bool +bernoulli p = (< p) <$> uniform +{-# INLINE bernoulli #-} + +binomial :: PrimMonad m => Int -> Double -> Prob m Int +binomial n p = liftM (length . filter id) $ replicateM n (bernoulli p) +{-# INLINE binomial #-} +