diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,5 @@
+# Revision history for dirichlet
+
+## 0.1.0.0 -- YYYY-mm-dd
+
+* First version. Released on an unsuspecting world.
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2020, Dominik Schrempf
+
+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 Dominik Schrempf 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/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/dirichlet.cabal b/dirichlet.cabal
new file mode 100644
--- /dev/null
+++ b/dirichlet.cabal
@@ -0,0 +1,43 @@
+cabal-version:       2.4
+
+name:                dirichlet
+version:             0.1.0.0
+synopsis:            Multivariate dirichlet distribution
+description:         Please see the README on GitHub at <https://github.com/dschrempf/dirichlet#readme>
+homepage:            https://github.com/dschrempf/dirichlet
+bug-reports:
+license:             BSD-3-Clause
+license-file:        LICENSE
+author:              Dominik Schrempf
+maintainer:          dominik.schrempf@gmail.com
+-- copyright:
+category:            Math
+build-type:          Simple
+extra-source-files:  CHANGELOG.md
+
+library
+  exposed-modules:     Statistics.Distribution.Dirichlet
+  -- other-modules:
+  -- other-extensions:
+  build-depends:       base >=4.14 && <4.15
+                     , log-domain
+                     , math-functions
+                     , mwc-random
+                     , primitive
+                     , vector
+  hs-source-dirs:      src
+  default-language:    Haskell2010
+
+test-suite dirichlet-test
+  default-language:    Haskell2010
+  type:                exitcode-stdio-1.0
+  hs-source-dirs:      test
+  main-is:             Spec.hs
+  other-modules:       Statistics.Distribution.DirichletSpec
+  build-depends:       base >=4.14 && <4.15
+                     , dirichlet
+                     , hspec
+                     , hspec-discover
+                     , log-domain
+                     , mwc-random
+                     , vector
diff --git a/src/Statistics/Distribution/Dirichlet.hs b/src/Statistics/Distribution/Dirichlet.hs
new file mode 100644
--- /dev/null
+++ b/src/Statistics/Distribution/Dirichlet.hs
@@ -0,0 +1,151 @@
+-- |
+-- Module      :  Statistics.Distribution.Dirichlet
+-- Description :  Multivariate Dirichlet distribution
+-- Copyright   :  (c) Dominik Schrempf, 2020
+-- License     :  GPL-3.0-or-later
+--
+-- Maintainer  :  dominik.schrempf@gmail.com
+-- Stability   :  unstable
+-- Portability :  portable
+--
+-- Creation date: Tue Oct 20 10:10:39 2020.
+module Statistics.Distribution.Dirichlet
+  ( -- * Dirichlet distribution
+    DirichletDistribution (ddGetParameters),
+    dirichletDistribution,
+    dirichletDensity,
+    dirichletSample,
+
+    -- * Symmetric Dirichlet distribution
+    DirichletDistributionSymmetric (ddSymGetParameter),
+    dirichletDistributionSymmetric,
+    dirichletDensitySymmetric,
+    dirichletSampleSymmetric,
+  )
+where
+
+import Control.Monad.Primitive
+import qualified Data.Vector.Unboxed as V
+import Numeric.Log
+import Numeric.SpecFunctions
+import System.Random.MWC
+import System.Random.MWC.Distributions
+
+-- | The Dirichlet distribution is identified by a vector of parameter values.
+data DirichletDistribution = DirichletDistribution
+  { ddGetParameters :: V.Vector Double,
+    _getDimension :: Int,
+    _getNormConst :: Log Double
+  }
+  deriving (Eq, Show)
+
+-- Check if vector is strictly positive.
+isNegativeOrZero :: V.Vector Double -> Bool
+isNegativeOrZero = V.any (<= 0)
+
+-- Inverse multivariate beta function. Does not check if parameters are valid!
+invBeta :: V.Vector Double -> Log Double
+invBeta v = Exp $ logDenominator - logNominator
+  where
+    logNominator = V.sum $ V.map logGamma v
+    logDenominator = logGamma (V.sum v)
+
+-- | Create a Dirichlet distribution from the given parameter vector.
+--
+-- Return Left if:
+-- - The parameter vector has less then two elements.
+-- - One or more parameters are negative or zero.
+dirichletDistribution :: V.Vector Double -> Either String DirichletDistribution
+dirichletDistribution v
+  | V.length v < 2 =
+    Left "dirichletDistribution: Parameter vector is too short."
+  | isNegativeOrZero v =
+    Left "dirichletDistribution: One or more parameters are negative or zero."
+  | otherwise = Right $ DirichletDistribution v (V.length v) (invBeta v)
+
+-- Tolerance.
+eps :: Double
+eps = 1e-14
+
+-- Check if vector is normalized with tolerance 'eps'.
+isNormalized :: V.Vector Double -> Bool
+isNormalized v
+  | abs (V.sum v - 1.0) > eps = False
+  | otherwise = True
+
+-- | Density of the Dirichlet distribution evaluated at a given value vector.
+--
+-- Return 0 if:
+-- - The value vector has a different length than the parameter vector.
+-- - The value vector has elements being negative or zero.
+-- - The value vector does not sum to 1.0 (with tolerance @eps = 1e-14@).
+dirichletDensity :: DirichletDistribution -> V.Vector Double -> Log Double
+dirichletDensity (DirichletDistribution as k c) xs
+  | k /= V.length xs = 0
+  | isNegativeOrZero xs = 0
+  | not (isNormalized xs) = 0
+  | otherwise = c * Exp logXsPow
+  where
+    logXsPow = V.sum $ V.zipWith (\a x -> log $ x ** (a - 1.0)) as xs
+
+-- | Sample a value vector from the Dirichlet distribution.
+dirichletSample :: PrimMonad m => DirichletDistribution -> Gen (PrimState m) -> m (V.Vector Double)
+dirichletSample (DirichletDistribution as _ _) g = do
+  ys <- V.mapM (\a -> gamma a 1.0 g) as
+  let s = V.sum ys
+  return $ V.map (/ s) ys
+
+-- | The Dirichlet distribution is identified by a vector of parameter values.
+data DirichletDistributionSymmetric = DirichletDistributionSymmetric
+  { ddSymGetParameter :: Double,
+    _symGetDimension :: Int,
+    _symGetNormConst :: Log Double
+  }
+  deriving (Eq, Show)
+
+-- Inverse multivariate beta function. Does not check if parameters are valid!
+invBetaSym :: Int -> Double -> Log Double
+invBetaSym k a = Exp $ logDenominator - logNominator
+  where
+    logNominator = fromIntegral k * logGamma a
+    logDenominator = logGamma (fromIntegral k * a)
+
+-- | Create a symmetric Dirichlet distribution of given dimension and parameter.
+--
+-- Return Left if:
+-- - The given dimension is smaller than two.
+-- - The parameter is negative or zero.
+dirichletDistributionSymmetric :: Int -> Double -> Either String DirichletDistributionSymmetric
+dirichletDistributionSymmetric k a
+  | k < 2 =
+    Left "dirichletDistributionSymmetric: The dimension is smaller than two."
+  | a <= 0 =
+    Left "dirichletDistributionSymmetric: The parameter is negative or zero."
+  | otherwise = Right $ DirichletDistributionSymmetric a k (invBetaSym k a)
+
+-- | Density of the symmetric Dirichlet distribution evaluated at a given value
+-- vector.
+--
+-- Return 0 if:
+-- - The value vector has a different dimension.
+-- - The value vector has elements being negative or zero.
+-- - The value vector does not sum to 1.0 (with tolerance @eps = 1e-14@).
+dirichletDensitySymmetric :: DirichletDistributionSymmetric -> V.Vector Double -> Log Double
+dirichletDensitySymmetric (DirichletDistributionSymmetric a k c) xs
+  | k /= V.length xs = 0
+  | isNegativeOrZero xs = 0
+  | not (isNormalized xs) = 0
+  | otherwise = c * Exp logXsPow
+  where
+    logXsPow = V.sum $ V.map (\x -> log $ x ** (a - 1.0)) xs
+
+-- | Sample a value vector from the symmetric Dirichlet distribution.
+dirichletSampleSymmetric ::
+  PrimMonad m =>
+  DirichletDistributionSymmetric ->
+  Gen (PrimState m) ->
+  m (V.Vector Double)
+dirichletSampleSymmetric (DirichletDistributionSymmetric a k _) g = do
+  ys <- V.replicateM k (gamma a 1.0 g)
+  let s = V.sum ys
+  return $ V.map (/ s) ys
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 #-}
diff --git a/test/Statistics/Distribution/DirichletSpec.hs b/test/Statistics/Distribution/DirichletSpec.hs
new file mode 100644
--- /dev/null
+++ b/test/Statistics/Distribution/DirichletSpec.hs
@@ -0,0 +1,80 @@
+-- |
+-- Module      :  DirichletSpec
+-- Description :  Unit tests for DirichletSpec
+-- Copyright   :  (c) Dominik Schrempf, 2020
+-- License     :  GPL-3.0-or-later
+--
+-- Maintainer  :  dominik.schrempf@gmail.com
+-- Stability   :  unstable
+-- Portability :  portable
+--
+-- Creation date: Tue Oct 20 10:40:06 2020.
+module Statistics.Distribution.DirichletSpec
+  ( spec,
+  )
+where
+
+import Control.Monad
+import Data.Either
+import qualified Data.Vector.Unboxed as V
+import Numeric.Log hiding (sum)
+import Statistics.Distribution.Dirichlet
+import Test.Hspec
+import System.Random.MWC
+
+eps :: Double
+eps = 1e-14
+
+dd3 :: DirichletDistribution
+dd3 = either error id $ dirichletDistribution $ V.fromList [0.2, 10, 20]
+
+alphas10 :: V.Vector Double
+alphas10 = V.fromList [0.1, 0.2, 0.3, 0.4, 0.5, 16, 17, 18, 19, 20]
+
+dd10 :: DirichletDistribution
+dd10 = either error id $ dirichletDistribution alphas10
+
+ddSym :: Int -> Double -> DirichletDistribution
+ddSym n a = either error id $ dirichletDistribution $ V.replicate n a
+
+-- Extract means.
+xbar :: Int -> [V.Vector Double] -> Double
+xbar n xss = sum xs / fromIntegral (length xs)
+  where xs = map (V.! n) xss
+
+spec :: Spec
+spec = do
+  describe "dirichletDistribution" $ do
+    it "only works for valid parameter vectors" $ do
+      dirichletDistribution (V.fromList []) `shouldSatisfy` isLeft
+      dirichletDistribution (V.fromList [1]) `shouldSatisfy` isLeft
+      dirichletDistribution (V.fromList [-3, 3, 192]) `shouldSatisfy` isLeft
+      dirichletDistribution (V.fromList [0, 3, 192]) `shouldSatisfy` isLeft
+  describe "dirichletDensity" $ do
+    it "is correct for some test cases" $ do
+      let ddSym2 = ddSym 2 0.5
+      let r = dirichletDensity ddSym2 (V.fromList [0.2, 0.8])
+      abs (exp (ln r) - 0.7957747154594766) `shouldSatisfy` (< eps)
+      let rBound = dirichletDensity ddSym2 (V.fromList [0, 1])
+      rBound `shouldBe` 0
+      let rOutBound1 = dirichletDensity ddSym2 (V.fromList [0, 1.1])
+      rOutBound1 `shouldBe` 0
+      let rOutBound2 = dirichletDensity ddSym2 (V.fromList [-0.1, 0.9])
+      rOutBound2 `shouldBe` 0
+      let r3 = dirichletDensity dd3 (V.fromList [0.3, 0.3, 0.4])
+      abs (exp (ln r3) - 0.0001217825570884453) `shouldSatisfy` (< eps)
+      let rWrongDim = dirichletDensity dd3 (V.fromList [0.3, 0.7])
+      rWrongDim `shouldBe` 0
+  describe "dirichletSample" $ do
+    it "returns valid value vectors with expected mean" $ do
+      g <- create
+      let ddSym10 = ddSym 10 10
+      xs <- replicateM 1000 (dirichletSample ddSym10 g)
+      map V.length xs `shouldBe` replicate 1000 10
+      -- print [ xbar i xs | i <- [0..9]]
+      [ abs (xbar i xs - 0.1) > 0.01 | i <- [0..9]] `shouldBe` replicate 10 False
+      xs <- replicateM 1000 (dirichletSample dd10 g)
+      map V.length xs `shouldBe` replicate 1000 10
+      -- print [ xbar i xs | i <- [0..9]]
+      let aSum = V.sum alphas10
+      [ abs (xbar i xs - (alphas10 V.! i / aSum)) > 0.01 | i <- [0..9]] `shouldBe` replicate 10 False