stochastic 0.1.0.0 → 0.1.1.0
raw patch · 15 files changed
+250/−91 lines, 15 filesdep −math-functionsdep ~basedep ~containersdep ~mtl
Dependencies removed: math-functions
Dependency ranges changed: base, containers, mtl, random
Files
- ChangeLog.md +6/−0
- src/Data/Stochastic.hs +70/−27
- src/Data/Stochastic/Internal.hs +28/−9
- src/Data/Stochastic/Types.hs +62/−15
- stochastic.cabal +13/−7
- test/Beta.hs +36/−0
- test/ChartTest.hs +4/−4
- test/ContrivedGambler.hs +4/−4
- test/CoolCharts.hs +6/−6
- test/MonadLaws.hs +3/−3
- test/MontyHall.hs +6/−4
- test/Normal10.hs +3/−3
- test/Normal3.hs +3/−3
- test/Swindler.hs +3/−3
- test/Test.hs +3/−3
ChangeLog.md view
@@ -3,3 +3,9 @@ ## 0.1.0.0 -- 2016-08-19 * First version. Released on an unsuspecting world.++## 0.1.1.0 -- 2016-08-26++* Exposes composeProcess for quickly composing StochProcesses.+* Un-reinvents the State monad.+* Adds gamma and beta distributions.
src/Data/Stochastic.hs view
@@ -18,11 +18,15 @@ module Data.Stochastic ( -- * Constructing a Sample- certain-, uniform+ mkSample+, certain+, discreteUniform , discrete , bernoulli , normal+, uniform+, gamma+, beta -- * Sampling from a Sample , sample , sample_@@ -31,10 +35,15 @@ , sampleIO_ , sampleION -- * Constructing a StochProcess+, toProcess , certainProcess-, uniformProcess+, discreteUniformProcess , discreteProcess , normalProcess+, uniformProcess+, gammaProcess+, betaProcess+, composeProcess -- * Running a StochProcess , runProcess , runProcess_@@ -51,6 +60,7 @@ , module Data.Stochastic.Internal ) where +import Control.Monad.State import Control.Monad.Trans import Control.Monad.Writer @@ -124,32 +134,42 @@ -- | 'StochProcess' sample for a normal distribution that records -- the value sampled from the normal distribution. normalProcess :: Mean -> StDev -> StochProcess -normalProcess mean std = do- sample <- lift $ normal mean std- tell $ S.singleton sample- return sample+normalProcess mean std = toProcess $ normal mean std --- | 'StochProcess' sample for a distribution over 'Double's that always+-- | 'StochProcess' for a distribution over 'Double's that always -- returns the same value when sampled, and records that value. certainProcess :: Double -> StochProcess -certainProcess a = do- sample <- lift $ certain a- tell $ S.singleton sample- return sample+certainProcess a = toProcess $ certain a --- | 'StochProcess' sample for a discrete distribution over 'Double's+-- | 'StochProcess' for a discrete distribution over 'Double's -- that records the value sampled from the normal distribution. discreteProcess :: [(Double, Double)] -> StochProcess -discreteProcess a = do- sample <- lift $ discrete a- tell $ S.singleton sample- return sample+discreteProcess a = toProcess $ discrete a --- | 'StochProcess' sample for a uniform distribution over 'Double's--- that records the value sampled from it.-uniformProcess :: [Double] -> StochProcess-uniformProcess l = do- sample <- lift $ uniform l+-- | 'StochProcess' for a discrete uniform distribution +-- over 'Double's that records the value sampled from it.+discreteUniformProcess :: [Double] -> StochProcess+discreteUniformProcess l = toProcess $ discreteUniform l++-- | 'StochProcess' for a discrete uniform distribution +-- over 'Double's that records the value sampled from it.+uniformProcess :: StochProcess+uniformProcess = toProcess uniform++-- | 'StochProcess' for a gamma distribution with+-- provided shape and scale parameters.+gammaProcess :: Double -> Double -> StochProcess+gammaProcess a b = toProcess $ gamma a b++-- | 'StochProcess' for a beta distribution.+betaProcess :: Double -> Double -> StochProcess+betaProcess a b = toProcess $ beta a b++-- | Function to create a 'StochProcess' out of a provided+-- 'Sample' over 'Double's.+toProcess :: Sample StdGen Distribution Double -> StochProcess+toProcess s = do+ sample <- lift s tell $ S.singleton sample return sample @@ -171,16 +191,39 @@ discrete [] = error "do not construct empty discrete distributions" discrete l = mkSample $ Discrete l --- | 'Sample' for a uniform distribution+-- | 'Sample' for a discrete uniform distribution -- given a list of provided values.-uniform :: (RandomGen g) => [a] -> Sample g Distribution a-uniform l = mkSample $ Uniform l+discreteUniform :: (RandomGen g) => [a] -> Sample g Distribution a+discreteUniform l = mkSample $ DiscreteUniform l +-- | 'Sample' for a continuous uniform distribution+-- with support [0, 1].+uniform :: (RandomGen g) => Sample g Distribution Double+uniform = mkSample Uniform+ -- | 'Sample' for a distribution where we always sample -- the same value. certain :: (RandomGen g, Sampleable d) => a -> Sample g d a certain = mkSample . certainDist +-- | 'Sample' for a gamma distribution given shape parameter+-- and scale parameter.+gamma :: RandomGen g+ => Double + -- ^ The shape parameter.+ -> Double + -- ^ The scale parameter.+ -> Sample g Distribution Double+gamma a b = if a <= 0 || b <= 0 + then error "cannot construct gamma dist with <=0 params"+ else mkSample $ Gamma a b++-- | 'Sample' for a beta distribution.+beta :: RandomGen g => Double -> Double -> Sample g Distribution Double+beta a b = if a <= 0 || b <= 0 + then error "cannot construct gamma dist with <=0 params"+ else mkSample $ Beta a b+ -- | Get one sample of type a from the 'Sample' along with -- a new 'StdGen'. --@@ -188,7 +231,7 @@ -- 'RandomGen' because when we sample from normal -- distributions, we consume one extra 'RandomGen'. sample :: (RandomGen g, Sampleable d) => Sample g d a -> g -> (a, g)-sample s g = let (dist, g') = runSample s g+sample s g = let (dist, g') = runState (runSample s) g (a, g'') = sampleFrom dist g' in (a, snd $ next g'') @@ -223,4 +266,4 @@ -- | Function to make a 'Sample' out of a provided -- 'Distribution'. mkSample :: (RandomGen g, Sampleable d) => d a -> Sample g d a-mkSample d = Sample $ \g -> (d, snd $ next g)+mkSample d = Sample $ return d
src/Data/Stochastic/Internal.hs view
@@ -8,21 +8,40 @@ -} module Data.Stochastic.Internal ( boxMuller-, decentRandom+, closedRnd+, openRnd+, closedOpenRnd+, openClosedRnd ) where import System.Random -import Numeric.MathFunctions.Constants- -- | Function to convert values sampled from uniform distribution -- to a value sampled from a standard normal distribution. boxMuller :: (Floating a) => a -> a -> a boxMuller u1 u2 = cos (2 * pi * u2) * (sqrt $ (-2) * (log u1)) --- | Randoms that aren't too small or equal to 1.-decentRandom :: (RandomGen g) => g -> (Double, g)-decentRandom gen = let (sampled, newG) = randomR (0, 1.0) gen- in if sampled <= m_epsilon || sampled == 1- then decentRandom newG- else (sampled, newG)+-- | Randoms in the interval [0, 1]+closedRnd :: (RandomGen g) => g -> (Double, g)+closedRnd gen = randomR (0, 1.0) gen++-- | Randoms in the interval (0, 1)+openRnd :: (RandomGen g) => g -> (Double, g)+openRnd gen = let (a, g) = closedRnd gen+ in if a == 0 || a == 1+ then openRnd g+ else (a, g)+ +-- | Randoms in the interval [0, 1)+closedOpenRnd :: (RandomGen g) => g -> (Double, g)+closedOpenRnd gen = let (a, g) = closedRnd gen+ in if a == 1+ then closedOpenRnd g+ else (a, g)+ +-- | Randoms in the interval (0, 1]+openClosedRnd :: (RandomGen g) => g -> (Double, g)+openClosedRnd gen = let (a, g) = closedRnd gen+ in if a == 0+ then openClosedRnd g+ else (a, g)
src/Data/Stochastic/Types.hs view
@@ -27,9 +27,11 @@ , Sampler (..) , Mean (..) , StDev (..)+, marsagliaTsang ) where import Control.Monad+import Control.Monad.State import Control.Monad.Writer import Data.Stochastic.Internal@@ -46,8 +48,14 @@ Normal :: Mean -> StDev -> Distribution Double Bernoulli :: Double -> Distribution Bool Discrete :: [(a, Double)] -> Distribution a- Uniform :: [a] -> Distribution a+ DiscreteUniform :: [a] -> Distribution a+ Uniform :: Distribution Double Certain :: a -> Distribution a+ Gamma :: Double -> Double -> Distribution Double+ -- ^ Gamma distribution, where the first parameter is the+ -- shape parameter alpha, and the second parameter is the+ -- scale parameter beta.+ Beta :: Double -> Double -> Distribution Double -- | Class of types from which samples can be obtained. class Sampleable d where@@ -65,15 +73,17 @@ sampleFrom da g = case da of Normal mean stdev - -> let (a, g') = decentRandom g - (a', g'') = decentRandom g'+ -> let (a, g') = closedRnd g + (a', g'') = closedRnd g' s = (stdev * (boxMuller a a')) + mean in (s, g') Bernoulli prob - -> let (a, g') = decentRandom g+ -> let (a, g') = closedRnd g in (a <= prob, g')+ Discrete []+ -> error "cannot sample from empty discrete distribution" Discrete l - -> let (a, g') = decentRandom g+ -> let (a, g') = closedRnd g in (scan a l, g') where scan lim [] = if lim <= 0 then error $ "not normalized discrete dist"@@ -81,28 +91,60 @@ scan lim (x:xs) = if lim <= snd x then fst x else scan (lim - snd x) xs- Uniform l- -> let (a, g') = decentRandom g+ DiscreteUniform []+ -> error "cannot sample from empty discrete distribution"+ DiscreteUniform l+ -> let (a, g') = closedRnd g prob = 1 / (fromIntegral $ length l) in (l !! (floor $ a / prob), g')+ Uniform+ -> closedRnd g+ Gamma alpha beta+ -> if alpha <= 0 || beta <= 0 then error "alpha and beta parameter cannot be less than or equal to zero in beta distribution"+ else if alpha > 0 && alpha < 1 then+ let (a, g') = sampleFrom (Gamma (alpha + 1) beta) g+ (uni, g'') = openRnd g'+ in (a * (uni ** (1/alpha)), g'')+ else let d = alpha - (1/3)+ c = 1 / sqrt (9 * d)+ (m, g') = marsagliaTsang d c g+ in (m * beta, g')+ Beta alpha beta+ -> let (x, g') = sampleFrom (Gamma alpha 1) g+ (y, g'') = sampleFrom (Gamma beta 1) g'+ in (x / (x + y), g'') Certain val - -> (val, snd $ decentRandom g) + -> (val, snd $ openRnd g) -- Seemingly unnecessary, but important to obey the monad laws to always produce the same RandomGen each time we sample. certainDist = Certain +-- | Marsaglia and Tsang's rejection method+-- for generating Gamma variates with parameters+-- alpha and 1, where 1 is the scale parameter,+-- given d and c.+marsagliaTsang :: (RandomGen g) => Double -> Double -> g -> (Double, g)+marsagliaTsang d c g =+ let (norm, g') = sampleFrom (Normal 0 1) g+ (uni, g'') = sampleFrom Uniform g'+ v = (1 + (c * norm)) ** 3+ in if norm > ((-1)/c) && + log uni < ((norm ** 2)/2 + d - (d * v) + (d * log v))+ then (d * v, g'')+ else marsagliaTsang d c g''+ -- | Show instance for 'Distribution's. instance (Show a) => Show (Distribution a) where show da = case da of Normal mean stdev -> "Normal " ++ show mean ++ " " ++ show stdev Bernoulli prob -> "Bernoulli " ++ show prob Discrete l -> "Discrete " ++ show l- Uniform l -> "Uniform " ++ show l + DiscreteUniform l -> "DiscreteUniform " ++ show l Certain val -> "Certain " ++ show val -- | 'Sample' monad containing a random number generator plus a type from which -- we can sample values of type a newtype Sample g d a- = Sample { runSample :: (RandomGen g, Sampleable d) => g -> (d a, g) }+ = Sample { runSample :: (RandomGen g, Sampleable d) => State g (d a) } -- | Monad that represents a stochastic process. -- It allows us to record numeric values as we sample.@@ -111,11 +153,16 @@ -- | Monad instance for Sample. instance (RandomGen g, Sampleable d) => Monad (Sample g d) where- return x = Sample $ \g -> (certainDist x, snd $ next g)- (>>=) ma f = Sample $ \g -> - let (dist, g') = runSample ma g- (a, g'') = sampleFrom dist g'- in runSample (f a) g''+ return x = Sample $ do+ modify (snd . next)+ return $ certainDist x++ (>>=) ma f = Sample $ do+ modify (snd . next)+ dist <- runSample ma+ g <- get+ let a = fst $ sampleFrom dist g+ runSample (f a) -- | Trivial 'Functor' instance for 'Sample' 'StdGen' 'Distribution'. instance (RandomGen g, Sampleable s) => Functor (Sample g s) where
stochastic.cabal view
@@ -10,7 +10,7 @@ -- PVP summary: +-+------- breaking API changes -- | | +----- non-breaking API additions -- | | | +--- code changes with no API change-version: 0.1.0.0+version: 0.1.1.0 -- A short (one-line) description of the package. synopsis: Monadic composition of probabilistic functions and sampling.@@ -22,7 +22,7 @@ can then be sampled from to generate datapoints. -- URL for the project homepage or repository.-homepage: http://kevinl.io/posts/2016-08-17-sampling-monad.html+homepage: http://kevinl.io/posts/2016-08-24-sampling-monad.html -- The license under which the package is released. license: GPL-3@@ -55,7 +55,7 @@ library -- Modules exported by the library. exposed-modules: Data.Stochastic- -- Data.Sample.Chart+ -- Data.Stochastic.Chart Data.Stochastic.Types Data.Stochastic.Internal @@ -71,7 +71,6 @@ -- Chart-cairo >=1.8 && <1.9, containers >= 0.5 && <0.6, random >=1.1 && <2,- math-functions >=0.2 && <0.3, mtl >=2.2 && <2.3 -- Directories containing source files.@@ -118,7 +117,7 @@ -- | Options to pass to ghc ghc-options: -O2 -threaded- "-with-rtsopts=-N -pa -s -h -i0.1"+ -- "-with-rtsopts=-N -pa -s -h -i0.1" test-suite normal10 type: exitcode-stdio-1.0 hs-source-dirs: test@@ -128,7 +127,7 @@ -- | Options to pass to ghc ghc-options: -O2 -threaded- "-with-rtsopts=-N -pa -s -h -i0.1"+ -- "-with-rtsopts=-N -pa -s -h -i0.1" test-suite chart type: exitcode-stdio-1.0@@ -163,7 +162,7 @@ type: exitcode-stdio-1.0 hs-source-dirs: test main-is: MontyHall.hs- build-depends: base, stochastic, containers, random+ build-depends: base, stochastic, containers, random, mtl default-language: Haskell2010 -- | Options to pass to ghc ghc-options: -O2@@ -178,3 +177,10 @@ -- | Options to pass to ghc ghc-options: -O2 -threaded++test-suite beta+ type: exitcode-stdio-1.0+ hs-source-dirs: test+ main-is: Beta.hs+ build-depends: base, stochastic, random, Chart-cairo, Chart, containers, mtl+ default-language: Haskell2010
+ test/Beta.hs view
@@ -0,0 +1,36 @@+import Control.Monad.Trans++import Data.Foldable+import Data.Stochastic+import Data.Stochastic.Types+import Data.Stochastic.Internal+import Data.Stochastic.Chart++import qualified Data.Sequence as S++import System.Random+import Graphics.Rendering.Chart.Backend.Cairo+import Graphics.Rendering.Chart.Easy++beta1 :: Sampler Double+beta1 = beta 0.5 0.5++beta2 :: Sampler Double+beta2 = beta 0.5 2++gamma1 :: Sampler Double+gamma1 = gamma 3 2++gamma2 :: Sampler Double+gamma2 = gamma 3 0.5++main = do+ gen <- newStdGen+ let betas1 = sampleN 100000 beta1 gen+ let betas2 = sampleN 100000 beta2 gen+ let gammas1 = sampleN 100000 gamma1 gen+ let gammas2 = sampleN 100000 gamma2 gen+ renderableToFile def "charts/beta1.png" $ toRenderable $ histogram "Beta1" $ toList betas1+ renderableToFile def "charts/beta2.png" $ toRenderable $ histogram "Beta2" $ toList betas2+ renderableToFile def "charts/gamma1.png" $ toRenderable $ histogram "Gamma1" $ toList gammas1+ renderableToFile def "charts/gamma2.png" $ toRenderable $ histogram "Gamma2" $ toList gammas2
test/ChartTest.hs view
@@ -2,10 +2,10 @@ import Control.Monad.Trans -import Data.Sample-import Data.Sample.Types-import Data.Sample.Lib-import Data.Sample.Chart+import Data.Stochastic+import Data.Stochastic.Types+import Data.Stochastic.Internal+import Data.Stochastic.Chart import qualified Data.Sequence as S
test/ContrivedGambler.hs view
@@ -1,9 +1,9 @@ import Data.Foldable -import Data.Sample-import Data.Sample.Lib-import Data.Sample.Types-import Data.Sample.Chart+import Data.Stochastic+import Data.Stochastic.Internal+import Data.Stochastic.Types+import Data.Stochastic.Chart import qualified Data.Sequence as S
test/CoolCharts.hs view
@@ -1,9 +1,9 @@ import Control.Monad.Trans -import Data.Sample-import Data.Sample.Types-import Data.Sample.Lib-import Data.Sample.Chart+import Data.Stochastic+import Data.Stochastic.Types+import Data.Stochastic.Internal+import Data.Stochastic.Chart import qualified Data.Sequence as S @@ -30,5 +30,5 @@ main = do gen <- newStdGen- let a = runProcess_ normal1 gen- toFile def "chart.png" $ testChart $ S.singleton a+ let a = runProcessN 10 normal1 gen+ toFile def "chart.png" $ testChart $ a
test/MonadLaws.hs view
@@ -1,6 +1,6 @@-import Data.Sample-import Data.Sample.Lib-import Data.Sample.Types+import Data.Stochastic+import Data.Stochastic.Internal+import Data.Stochastic.Types import System.Random import System.Exit
test/MontyHall.hs view
@@ -1,7 +1,9 @@-import Data.Sample-import Data.Sample.Lib-import Data.Sample.Types+import Control.Monad.State +import Data.Stochastic+import Data.Stochastic.Internal+import Data.Stochastic.Types+ import qualified Data.Sequence as S import System.Random@@ -13,7 +15,7 @@ doorInit n = uniform $ Car : (take (n-1) $ repeat Goat) removeDoor :: Prize -> Sampler Prize -> Sampler Prize-removeDoor b sa = let (da, g) = runSample sa $ mkStdGen 0+removeDoor b sa = let (da, g) = runState (runSample sa) (mkStdGen 0) in case da of Uniform l -> uniform $ removeFst b l _ -> error "not monty hall"
test/Normal10.hs view
@@ -1,6 +1,6 @@-import Data.Sample-import Data.Sample.Types-import Data.Sample.Lib+import Data.Stochastic+import Data.Stochastic.Types+import Data.Stochastic.Internal import System.Random
test/Normal3.hs view
@@ -1,6 +1,6 @@-import Data.Sample-import Data.Sample.Types-import Data.Sample.Lib+import Data.Stochastic+import Data.Stochastic.Types+import Data.Stochastic.Internal import System.Random
test/Swindler.hs view
@@ -1,6 +1,6 @@-import Data.Sample-import Data.Sample.Lib-import Data.Sample.Types+import Data.Stochastic+import Data.Stochastic.Internal+import Data.Stochastic.Types import qualified Data.Sequence as S
test/Test.hs view
@@ -1,8 +1,8 @@ module Main where -import Data.Sample-import Data.Sample.Lib-import Data.Sample.Types+import Data.Stochastic+import Data.Stochastic.Internal+import Data.Stochastic.Types import qualified Data.Sequence as S