lazyppl 1.0 → 1.0.1

raw patch · 10 files changed

+170/−75 lines, 10 filesdep −monad-extrasdep ~basedep ~containersdep ~hmatrixPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies removed: monad-extras

Dependency ranges changed: base, containers, hmatrix, log-domain, math-functions, mtl, random, transformers

API changes (from Hackage documentation)

Files

• .DS_Store binary
• ChangeLog.md +13/−0
• lazyppl.cabal +21/−16
• src/.DS_Store binary
• src/LazyPPL.hs +10/−3
• src/LazyPPL/.DS_Store binary
• src/LazyPPL/Distributions.hs +94/−22
• src/LazyPPL/Distributions/.DS_Store binary
• src/LazyPPL/Distributions/IBP.hs +31/−29
• src/LazyPPL/Distributions/Memoization.hs +1/−5

binary file changed (6148 → absent bytes)

@@ -0,0 +1,13 @@+# Changelog for lazyppl++## 1.0.1++- Add support for GHC 9.8 and 9.10+- Drop upper bounds on dependencies (except `base`)+- Add log-density / log-PMF functions: `normalLogPdf`, `gammaLogPdf`, etc.+- Export `scoreLog` and `scoreProductLog` (already defined in 1.0, now visible)+- Remove `monad-extras` dependency++## 1.0++- Initial Hackage release

@@ -1,23 +1,32 @@ cabal-version: 1.12 --- This file has been generated from package.yaml by hpack version 0.35.1.+-- This file has been generated from package.yaml by hpack version 0.38.1. -- -- see: https://github.com/sol/hpack  name:           lazyppl-version:        1.0+version:        1.0.1 synopsis:       Lazy Probabilistic Programming Library-description:    LazyPPL is a Haskell library for Bayesian probabilistic programming. It supports lazy use of probability, which is useful for specifying non-parametric models, and we provide new Metropolis-Hastings algorithms to allow this. For illustrations, see <https://lazyppl-team.github.io/>. The Gaussian Process module uses hmatrix, which requires a separate lapack installation (or -dev or -devel).+description:    LazyPPL is a Haskell library for Bayesian probabilistic programming. It supports lazy use of probability, which is useful for specifying non-parametric models, and we provide new Metropolis-Hastings algorithms to allow this. For illustrations, see <https://lazyppl-team.github.io/>.+                .+                Note: the Gaussian process module (@LazyPPL.Distributions.GP@) uses @hmatrix@, which links against LAPACK. On Ubuntu install with @sudo apt-get install libgsl0-dev liblapack-dev@. category:       Statistics stability:      experimental homepage:       https://lazyppl-team.github.io/ bug-reports:    https://github.com/lazyppl-team/lazyppl/issues author:         LazyPPL team maintainer:     sam.staton@cs.ox.ac.uk-copyright:      2021-2024 LazyPPL team+copyright:      2021-2026 LazyPPL team license:        MIT license-file:   LICENSE build-type:     Simple+tested-with:+    GHC == 8.10.7+  , GHC == 9.6.7+  , GHC == 9.8.4+  , GHC == 9.10.1+extra-source-files:+    ChangeLog.md  source-repository head   type: git@@ -32,19 +41,15 @@       LazyPPL.Distributions.GP       LazyPPL.Distributions.IBP       LazyPPL.Distributions.Memoization-  other-modules:-      Paths_lazyppl   hs-source-dirs:       src   build-depends:-      base >=4.7 && <5-    , containers >=0.6.5.1 && <0.7-    , hmatrix >=0.20.2 && <0.21-    , log-domain >=0.13.2 && <0.14-    , math-functions >=0.3.4.2 && <0.3.5-    , monad-extras >=0.6.0 && <0.7-    , mtl >=2.2.2 && <2.4-    , random >=1.2.0 && <1.3-    , transformers >=0.5.6.2 && <0.7+      base >=4.7 && <9+    , containers >=0.6.5.1+    , hmatrix >=0.20.2+    , log-domain >=0.13.2+    , math-functions >=0.3.4.2+    , mtl >=2.2.2+    , random >=1.2.0+    , transformers >=0.5.6.2   default-language: Haskell2010-

binary file changed (6148 → absent bytes)

@@ -1,4 +1,6 @@ {-# LANGUAGE GeneralizedNewtypeDeriving, ScopedTypeVariables #-}+{-# OPTIONS_GHC -Wno-unrecognised-pragmas #-}+{-# HLINT ignore "Use second" #-}  {- | LazyPPL is a library for Bayesian probabilistic programming. It supports lazy use of probability, and we provide new Metropolis-Hastings simulation algorithms to allow this. Laziness appears to be a good paradigm for non-parametric statistics.  @@ -48,7 +50,7 @@       -- * Basic interface       --       -- | There are three building blocks for measures: `uniform` for probability measures; `sample` and `score` for unnormalized measures. Combined with the monad structure, these give all s-finite measures.-      uniform, sample, score,+      uniform, sample, score, scoreLog, scoreProductLog,       -- * Monte Carlo simulation       --       -- | The `Meas` type describes unnormalized measures. Monte Carlo simulation allows us to sample from an unnormalized measure. Our main Monte Carlo simulator is `mh`. @@ -61,7 +63,6 @@ import Data.Monoid import System.Random hiding (uniform) import Control.Monad-import Control.Monad.Extra import Control.Monad.State.Lazy (State, state , put, get, runState) import Numeric.Log @@ -188,7 +189,7 @@     When 1/@p@ is roughly the number of used sites, then this will be a bit like "single-site lightweight" MH.     If @p@ = 1 then this is "multi-site lightweight" MH. -    __Tip:__ if @m :: Prob a@ then use @map fst <$> (mh 1 $ sample m)@ to get a stream of samples from a probability distribution without conditioning. +    __Tip:__ if @m :: Prob a@ then use @fmap (map fst) (mh 1 $ sample m)@ to get a stream of samples from a probability distribution without conditioning.  --} {-- The algorithm is as follows: @@ -330,4 +331,10 @@ every n xs = case drop (n -1) xs of   (y : ys) -> y : every n ys   [] -> []++-- | Monadic iterate (inlined from monad-extras).+iterateM :: Monad m => (a -> m a) -> a -> m [a]+iterateM f x = do+    x' <- f x+    (x':) `liftM` iterateM f x' 

binary file changed (6148 → absent bytes)

@@ -3,16 +3,23 @@ Sometimes both a distribution (type @Prob a@) and pdf (type @a -> Double@) are given. Distributions are useful for sampling, densities are used for scoring.   For more distributions, see the Statistics.Distribution in the statistics package. -}+{-# OPTIONS_GHC -Wno-unrecognised-pragmas #-}+{-# HLINT ignore "Avoid lambda using `infix`" #-}   module LazyPPL.Distributions (-       -- * Continuous distributions-       normal,normalPdf,exponential,expPdf,gamma, beta, dirichlet, uniformbounded,-       -- * Discrete distributions-       bernoulli, uniformdiscrete, categorical, poisson, poissonPdf,-       -- * Streams-       iid)-       where+        -- * Continuous distributions+        normal,normalPdf,normalLogPdf,exponential,expPdf,exponentialLogPdf,+        gamma, gammaLogPdf, beta, betaLogPdf, dirichlet, dirichletLogPdf,+        uniformbounded, uniformBoundedLogPdf,+        logNormalLogPdf,+        -- * Discrete distributions+        bernoulli, bernoulliLogPmf, uniformdiscrete, categorical, categoricalLogPmf,+        poisson, poissonPdf, poissonLogPmf,+        binomialLogPmf,+        -- * Streams+        iid)+      where  import LazyPPL (Prob,uniform) import Data.List (findIndex)@@ -23,14 +30,13 @@   [Normal distribution](https://en.wikipedia.org/wiki/Normal_distribution) -} normal :: Double -- ^ mu, mean-       -> Double -- ^ sigma, standard deviation-       -> Prob Double-normal m s = do +        -> Double -- ^ sigma, standard deviation+        -> Prob Double+normal m s = do   x <- uniform   return $ (- invErfc (2 * x)) * (m_sqrt_2 * s) + m- normalPdf :: Double -> Double -> Double -> Double-normalPdf m s x = exp ((-(x - m) * (x -m) / (2 * s * s)) - log (m_sqrt_2_pi * s))+normalPdf m s x = exp (- ((x - m) * (x - m) / (2 * s * s)) - log (m_sqrt_2_pi * s))   {-|@@ -38,12 +44,11 @@ -} exponential :: Double -- ^ lambda, rate             -> Prob Double-exponential rate = do +exponential rate = do   x <- uniform   return $ - (log x / rate)- expPdf :: Double -> Double -> Double-expPdf rate x = exp (-rate*x) * rate+expPdf rate x = exp (- (rate * x)) * rate  {-|   [Gamma distribution](https://en.wikipedia.org/wiki/Gamma_distribution)@@ -59,8 +64,8 @@   [Beta distribution](https://en.wikipedia.org/wiki/Beta_distribution) -} beta :: Double -- ^ alpha-     -> Double -- ^ beta-     -> Prob Double+      -> Double -- ^ beta+      -> Prob Double beta a b = do   x <- uniform   return $ invIncompleteBeta a b x@@ -77,8 +82,8 @@   return $ fromIntegral n  poissonPdf :: Double -> Integer -> Double-poissonPdf rate n = let result = exp(-rate) * rate ^^ (fromIntegral n) / (factorial (fromIntegral n)) in -  if (isInfinite result) || (isNaN result) then exp (-rate + (fromIntegral n) * log rate - logGamma (fromIntegral (n+1))) else result+poissonPdf rate n = let result = exp(-rate) * rate ^^ fromIntegral n / factorial (fromIntegral n) in+  if isInfinite result || isNaN result then exp (-rate + fromIntegral n * log rate - logGamma (fromIntegral (n+1))) else result   {-|@@ -94,8 +99,8 @@  -- | [Continuous uniform distribution on a bounded interval](https://en.wikipedia.org/wiki/Continuous_uniform_distribution) uniformbounded :: Double -- ^ lower-               -> Double -- ^ upper-               -> Prob Double+                -> Double -- ^ upper+                -> Prob Double uniformbounded lower upper = do   x <- uniform   return $ (upper - lower) * x + lower@@ -121,7 +126,7 @@ {-| [Categorical distribution](https://www.google.com/search?client=safari&rls=en&q=categorical+distribution&ie=UTF-8&oe=UTF-8): Takes a list of k numbers that sum to 1,      and returns a random number between 0 and (k-1), weighted accordingly -} categorical :: [Double] -> Prob Int-categorical xs = do +categorical xs = do   r <- uniform   case findIndex (>r) $ tail $ scanl (+) 0 xs of     Just i -> return i@@ -130,3 +135,70 @@ {-| Returns an infinite stream of samples from the given distribution. --} iid :: Prob a -> Prob [a] iid p = do r <- p; rs <- iid p; return $ r : rs++normalLogPdf :: Double -> Double -> Double -> Double+normalLogPdf mu sigma x =+  let d = x - mu+  in -0.5 * log (2 * pi) - log sigma - (d * d) / (2 * sigma * sigma)++exponentialLogPdf :: Double -> Double -> Double+exponentialLogPdf rate x+  | x < 0    = -1e300+  | otherwise = log rate - rate * x++gammaLogPdf :: Double -> Double -> Double -> Double+gammaLogPdf k theta x+  | x <= 0   = -1e300+  | otherwise = (k - 1) * log x - x / theta - k * log theta - logGamma k++betaLogPdf :: Double -> Double -> Double -> Double+betaLogPdf a b x+  | x <= 0 || x >= 1 = -1e300+  | otherwise = (a - 1) * log x + (b - 1) * log (1 - x) - logBeta a b++bernoulliLogPmf :: Double -> Bool -> Double+bernoulliLogPmf p True  = log (max 1e-300 p)+bernoulliLogPmf p False = log (max 1e-300 (1 - p))++categoricalLogPmf :: [Double] -> Int -> Double+categoricalLogPmf probs k+  | k < 0 || k >= length probs = -1e300+  | otherwise = log (max 1e-300 (probs !! k))++uniformBoundedLogPdf :: Double -> Double -> Double -> Double+uniformBoundedLogPdf lower upper x+  | x < lower || x > upper = -1e300+  | otherwise = -log (upper - lower)++dirichletLogPdf :: [Double] -> [Double] -> Double+dirichletLogPdf alphas xs+  | length alphas /= length xs = -1e300+  | any (<= 0) xs              = -1e300+  | otherwise =+      let lnB = sum (map logGamma alphas) - logGamma (sum alphas)+      in sum (zipWith (\a x -> (a - 1) * log x) alphas xs) - lnB++logNormalLogPdf :: Double -> Double -> Double -> Double+logNormalLogPdf mu sigma x+  | x <= 0   = -1e300+  | otherwise =+      let lx = log x+          d  = lx - mu+      in -0.5 * log (2 * pi) - log sigma - lx - (d * d) / (2 * sigma * sigma)++binomialLogPmf :: Int -> Double -> Int -> Double+binomialLogPmf n p k+  | k < 0 || k > n  = -1e300+  | p <= 0 && k > 0 = -1e300+  | p >= 1 && k < n = -1e300+  | otherwise =+      logGamma (fromIntegral n + 1)+      - logGamma (fromIntegral k + 1)+      - logGamma (fromIntegral (n - k) + 1)+      + fromIntegral k * log (max 1e-300 p)+      + fromIntegral (n - k) * log (max 1e-300 (1 - p))++poissonLogPmf :: Double -> Integer -> Double+poissonLogPmf lambda k+  | k < 0    = -1e300+  | otherwise = fromIntegral k * log lambda - lambda - logGamma (fromIntegral k + 1)

binary file changed (6148 → absent bytes)

@@ -1,4 +1,6 @@ {-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# OPTIONS_GHC -Wno-unrecognised-pragmas #-}+{-# HLINT ignore "Use null" #-}  {-| An implementation of the Indian buffet process by [Griffiths and Ghahramani](https://papers.nips.cc/paper_files/paper/2005/file/2ef35a8b78b572a47f56846acbeef5d3-Paper.pdf).@@ -6,7 +8,7 @@ We are using abstract types to hide the implementation details, inspired by [Exchangeable Random Processes and Data Abstraction](https://www.cs.ox.ac.uk/people/hongseok.yang/paper/pps17a.pdf).   Illustration: [Feature extraction example](https://lazyppl-team.github.io/AdditiveClusteringDemo.html). --} +-}  module LazyPPL.Distributions.IBP where @@ -18,40 +20,40 @@ import Data.List  - + -- Some abstract types -newtype Restaurant = R ([[Bool]], Counter)  +newtype Restaurant = R ([[Bool]], Counter) newtype Dish = D Int  deriving (Eq,Ord,Show,MonadMemo Prob)  newCustomer :: Restaurant -> Prob [Dish]-newCustomer (R (matrix, ref)) = do -    i <- readAndIncrement ref +newCustomer (R (matrix, ref)) = do+    i <- readAndIncrement ref     return [ D k | k <- [0..(length (matrix!!i) - 1)], matrix!!i!!k ] -            -newRestaurant :: Double -> Prob Restaurant ++newRestaurant :: Double -> Prob Restaurant newRestaurant alpha = do-        r <- uniform +        r <- uniform         ref <- newCounter-        matrix <- ibp alpha  -        return $ R (matrix, ref) +        matrix <- ibp alpha+        return $ R (matrix, ref)  -matrix :: Double -> Int -> [Int] -> Prob [[Bool]] -matrix alpha index features = -     do +matrix :: Double -> Int -> [Int] -> Prob [[Bool]]+matrix alpha index features =+     do         let i = fromIntegral index-        existingDishes <- mapM (\m -> bernoulli ((fromIntegral m) / i)) features -        let newFeatures = zipWith (\a -> \b -> if b then a + 1 else a) features existingDishes -        nNewDishes     <- fmap fromIntegral $ poisson (alpha / i) -        let fixZero = if features == [] && nNewDishes == 0 then 1 else nNewDishes  -        let newRow = existingDishes ++ (take fixZero $ repeat True) -        rest           <- matrix alpha (index + 1) (newFeatures ++ (take fixZero $ repeat 1)) -        return $ newRow : rest  +        existingDishes <- mapM (\m -> bernoulli (fromIntegral m / i)) features+        let newFeatures = zipWith (\ a b -> if b then a + 1 else a) features existingDishes+        nNewDishes     <- fromIntegral <$> poisson (alpha / i)+        let fixZero = if features == [] && nNewDishes == 0 then 1 else nNewDishes+        let newRow = existingDishes ++ take fixZero (repeat True)+        rest           <- matrix alpha (index + 1) (newFeatures ++ take fixZero (repeat 1))+        return $ newRow : rest  -- the distribution on matrices -ibp :: Double -> Prob [[Bool]]  -ibp alpha = matrix alpha 1 [] +ibp :: Double -> Prob [[Bool]]+ibp alpha = matrix alpha 1 []   @@ -68,17 +70,17 @@ Daniel M. Roy, Vikash Mansinghka, Noah Goodman, and Joshua Tenenbaum ICML Workshop on Nonparametric Bayesian, 2008.  --}-data RestaurantS = RS [Double] +newtype RestaurantS = RS [Double] -data DishS = DS Int deriving (Eq,Ord,Show)+newtype DishS = DS Int deriving (Eq,Ord,Show)  newCustomerS :: RestaurantS -> Prob [DishS]-newCustomerS (RS rs) = -    do fs <- mapM bernoulli rs-       return $ map DS $ findIndices id fs+newCustomerS (RS rs) = do +  fs <- mapM bernoulli rs+  return $ map DS $ findIndices id fs -newRestaurantS :: Double -> Prob RestaurantS -newRestaurantS a = fmap RS $ stickScale 1+newRestaurantS :: Double -> Prob RestaurantS+newRestaurantS a = RS <$> stickScale 1   where stickScale p = do r' <- beta a 1                           let r = p * r'                           -- Truncate when the probabilities are getting small

@@ -25,13 +25,9 @@ module LazyPPL.Distributions.Memoization (MonadMemo, memoize, generalmemoize, memrec) where  import Control.Monad-import Control.Monad.Extra-import Control.Monad.State.Lazy (State, get, put, runState, state) import Data.IORef-import Data.List import Data.Map (empty, insert, keys, lookup, size)-import Debug.Trace-import LazyPPL +import LazyPPL import System.IO.Unsafe  {-| Type class for memoizable argument types @a@ under a monad @m@ -}