packages feed

dirichlet (empty) → 0.1.0.0

raw patch · 7 files changed

+312/−0 lines, 7 filesdep +basedep +dirichletdep +hspecsetup-changed

Dependencies added: base, dirichlet, hspec, hspec-discover, log-domain, math-functions, mwc-random, primitive, vector

Files

+ CHANGELOG.md view
@@ -0,0 +1,5 @@+# Revision history for dirichlet++## 0.1.0.0 -- YYYY-mm-dd++* First version. Released on an unsuspecting world.
+ LICENSE view
@@ -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.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ dirichlet.cabal view
@@ -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
+ src/Statistics/Distribution/Dirichlet.hs view
@@ -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
+ test/Spec.hs view
@@ -0,0 +1,1 @@+{-# OPTIONS_GHC -F -pgmF hspec-discover #-}
+ test/Statistics/Distribution/DirichletSpec.hs view
@@ -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