diff --git a/ChangeLog.md b/ChangeLog.md
new file mode 100644
--- /dev/null
+++ b/ChangeLog.md
@@ -0,0 +1,3 @@
+# Changelog for splitmix-distributions
+
+## Unreleased changes
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright Author name here (c) 2021
+
+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 Author name here 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,31 @@
+# splitmix-distributions
+
+Random samplers for some common distributions, as well as a convenient interface for composing them, based on `splitmix`.
+
+
+## Usage
+
+Compose your random sampler out of simpler ones thanks to the Applicative and Monad interface, e.g. this is how you would declare and sample a binary mixture of Gaussian random variables:
+
+
+    import Control.Monad (replicateM)
+    import System.Random.SplitMix.Distributions (Gen, sample, bernoulli, normal)
+
+    process :: Gen Double
+    process = do
+        coin <- bernoulli 0.7
+        if coin
+        then
+            normal 0 2
+        else
+            normal 3 1
+
+    dataset :: [Double]
+    dataset = sample 1234 $ replicateM 20 process
+
+
+and sample your data in a pure (`sample`) or monadic (`sampleT`) setting.
+
+## Implementation details
+
+The library is built on top of `splitmix`, so the caveats on safety and performance that apply there are relevant here as well.
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/splitmix-distributions.cabal b/splitmix-distributions.cabal
new file mode 100644
--- /dev/null
+++ b/splitmix-distributions.cabal
@@ -0,0 +1,51 @@
+name:           splitmix-distributions
+version:        0.1.0.0
+description:    Random samplers for some common distributions, as well as a convenient interface for composing them, based on splitmix. Please see the README on GitHub at <https://github.com/ocramz/splitmix-distributions#readme>
+homepage:       https://github.com/ocramz/splitmix-distributions#readme
+bug-reports:    https://github.com/ocramz/splitmix-distributions/issues
+category:       Math
+synopsis:       Random samplers for some common distributions, based on splitmix.
+author:         Marco Zocca
+maintainer:     ocramz
+copyright:      2021 Marco Zocca
+license:        BSD3
+license-file:   LICENSE
+build-type:     Simple
+cabal-version:  1.12
+tested-with:    GHC == 8.10.4
+extra-source-files:
+    README.md
+    ChangeLog.md
+
+source-repository head
+  type: git
+  location: https://github.com/ocramz/splitmix-distributions
+
+library
+  exposed-modules:
+      System.Random.SplitMix.Distributions
+  hs-source-dirs:
+      src
+  build-depends:
+      base >=4.7 && <5
+    , erf
+    , mtl
+    , splitmix
+    , transformers
+  default-language: Haskell2010
+
+test-suite splitmix-distributions-test
+  type: exitcode-stdio-1.0
+  main-is: Spec.hs
+  hs-source-dirs:
+      test
+  ghc-options: -threaded -rtsopts -with-rtsopts=-N
+  build-depends:
+      base >=4.7 && <5
+    , erf
+    , hspec
+    , mtl
+    , splitmix
+    , splitmix-distributions
+    , transformers
+  default-language: Haskell2010
diff --git a/src/System/Random/SplitMix/Distributions.hs b/src/System/Random/SplitMix/Distributions.hs
new file mode 100644
--- /dev/null
+++ b/src/System/Random/SplitMix/Distributions.hs
@@ -0,0 +1,200 @@
+{-# language GeneralizedNewtypeDeriving #-}
+{-# options_ghc -Wno-unused-imports #-}
+{-|
+Random samplers for few common distributions, with an interface similar to that of @mwc-probability@.
+
+= Usage
+
+Compose your random sampler out of simpler ones thanks to the Applicative and Monad interface, e.g. this is how you would declare and sample a binary mixture of Gaussian random variables:
+
+@
+import Control.Monad (replicateM)
+import System.Random.SplitMix.Distributions (Gen, sample, bernoulli, normal)
+
+process :: `Gen` Double
+process = do
+  coin <- `bernoulli` 0.7
+  if coin
+    then
+      `normal` 0 2
+    else
+      normal 3 1
+
+dataset :: [Double]
+dataset = `sample` 1234 $ replicateM 20 process
+@
+
+and sample your data in a pure (`sample`) or monadic (`sampleT`) setting.
+
+== Implementation details
+
+The library is built on top of @splitmix@, so the caveats on safety and performance that apply there are relevant here as well.
+
+
+-}
+module System.Random.SplitMix.Distributions (
+  -- * Distributions
+  -- ** Continuous
+  stdUniform, uniformR,
+  exponential,
+  stdNormal, normal,
+  beta,
+  gamma,
+  -- ** Discrete
+  bernoulli,
+  -- * PRNG
+  -- ** Pure
+  Gen, sample,
+  -- ** Monadic
+  GenT, sampleT,
+  withGen
+                                            ) where
+
+import Control.Monad (replicateM)
+import Control.Monad.IO.Class (MonadIO(..))
+import Data.Functor.Identity (Identity(..))
+import GHC.Word (Word64)
+
+-- erf
+import Data.Number.Erf (InvErf(..))
+-- mtl
+import Control.Monad.Trans.Class (MonadTrans(..))
+import Control.Monad.State (MonadState(..), modify)
+-- splitmix
+import System.Random.SplitMix (SMGen, mkSMGen, splitSMGen, nextInt, nextInteger, nextDouble)
+-- transformers
+import Control.Monad.Trans.State (StateT(..), runStateT, evalStateT, State, runState, evalState)
+
+-- | Random generator
+--
+-- wraps 'splitmix' state-passing inside a 'StateT' monad
+--
+-- useful for embedding random generation inside a larger effect stack
+newtype GenT m a = GenT { unGen :: StateT SMGen m a } deriving (Functor, Applicative, Monad, MonadState SMGen, MonadTrans, MonadIO)
+
+-- | Pure random generation
+type Gen = GenT Identity
+
+-- | Monadic evaluation
+sampleT :: Monad m =>
+            Word64 -- ^ random seed
+         -> GenT m a -> m a
+sampleT seed gg = evalStateT (unGen gg) (mkSMGen seed)
+
+-- | Pure evaluation
+sample :: Word64 -- ^ random seed
+        -> Gen a
+        -> a
+sample seed gg = evalState (unGen gg) (mkSMGen seed)
+
+
+-- | Bernoulli trial
+bernoulli :: Double -- ^ bias parameter \( 0 \lt p \lt 1 \)
+          -> Gen Bool
+bernoulli p = withGen (bernoulliF p)
+
+-- | Uniform between two values
+uniformR :: Double -- ^ low
+         -> Double -- ^ high
+         -> Gen Double
+uniformR lo hi = scale <$> stdUniform
+  where
+    scale x = x * (hi - lo) + lo
+
+-- | Standard normal
+stdNormal :: Gen Double
+stdNormal = normal 0 1
+
+-- | Uniform in [0, 1)
+stdUniform :: Gen Double
+stdUniform = withGen nextDouble
+
+-- | Beta distribution, from two standard uniform samples
+beta :: Double -- ^ shape parameter \( \alpha \gt 0 \) 
+     -> Double -- ^ shape parameter \( \beta \gt 0 \)
+     -> Gen Double
+beta a b = go
+  where
+    go = do
+      (y1, y2) <- sample2
+      if
+        y1 + y2 <= 1
+        then pure (y1 / (y1 + y2))
+        else go
+    sample2 = f <$> stdUniform <*> stdUniform
+      where
+        f u1 u2 = (u1 ** (1/a), u2 ** (1/b))
+
+-- | Gamma distribution, using Ahrens-Dieter accept-reject (algorithm GD):
+--
+-- Ahrens, J. H.; Dieter, U (January 1982). "Generating gamma variates by a modified rejection technique". Communications of the ACM. 25 (1): 47–54
+gamma :: Double -- ^ shape parameter \( k \gt 0 \)
+      -> Double -- ^ scale parameter \( \theta \gt 0 \)
+      -> Gen Double
+gamma k th = do
+  xi <- sampleXi
+  us <- replicateM n (log <$> stdUniform)
+  pure $ th * xi - sum us
+  where
+    sampleXi = do
+      (xi, eta) <- sample2
+      if eta > xi ** (delta - 1) * exp (- xi)
+        then sampleXi
+        else pure xi
+    (n, delta) = (floor k, k - fromIntegral n)
+    ee = exp 1
+    sample2 = f <$> stdUniform <*> stdUniform <*> stdUniform
+      where
+        f u v w
+          | u <= ee / (ee + delta) =
+            let xi = v ** (1/delta)
+            in (xi, w * xi ** (delta - 1))
+          | otherwise =
+            let xi = 1 - log v
+            in (xi, w * exp (- xi))
+
+
+-- | Normal distribution
+normal :: Double -- ^ mean
+       -> Double -- ^ standard deviation \( \sigma \gt 0 \)
+       -> Gen Double
+normal mu sig = withGen (normalF mu sig)
+
+-- | Exponential distribution
+exponential :: Double -- ^ rate parameter \( \lambda > 0 \)
+            -> Gen Double
+exponential l = withGen (exponentialF l)
+
+-- | Wrap a 'splitmix' PRNG function
+withGen :: Monad m =>
+           (SMGen -> (a, SMGen)) -- ^ explicit generator passing (e.g. 'nextDouble')
+        -> GenT m a
+withGen f = GenT $ do
+  gen <- get
+  let
+    (b, gen') = f gen
+  put gen'
+  pure b
+
+exponentialF :: Double -> SMGen -> (Double, SMGen)
+exponentialF l g = (exponentialICDF l x, g') where (x, g') = nextDouble g
+
+normalF :: Double -> Double -> SMGen -> (Double, SMGen)
+normalF mu sig g = (normalICDF mu sig x, g') where (x, g') = nextDouble g
+
+bernoulliF :: Double -> SMGen -> (Bool, SMGen)
+bernoulliF p g = (x < p , g') where (x, g') = nextDouble g
+
+
+-- | inverse CDF of normal rv
+normalICDF :: InvErf a =>
+              a -- ^ mean
+           -> a -- ^ std dev
+           -> a -> a
+normalICDF mu sig p = mu + sig * sqrt 2 * inverf (2 * p - 1)
+
+-- | inverse CDF of exponential rv
+exponentialICDF :: Floating a =>
+                   a -- ^ rate
+                -> a -> a
+exponentialICDF l p = (- 1 / l) * log (1 - p)
diff --git a/test/Spec.hs b/test/Spec.hs
new file mode 100644
--- /dev/null
+++ b/test/Spec.hs
@@ -0,0 +1,1 @@
+{-# OPTIONS_GHC -F -pgmF hspec-discover #-}