diff --git a/.DS_Store b/.DS_Store
deleted file mode 100644
Binary files a/.DS_Store and /dev/null differ
diff --git a/ChangeLog.md b/ChangeLog.md
new file mode 100644
--- /dev/null
+++ b/ChangeLog.md
@@ -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
diff --git a/lazyppl.cabal b/lazyppl.cabal
--- a/lazyppl.cabal
+++ b/lazyppl.cabal
@@ -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
-
diff --git a/src/.DS_Store b/src/.DS_Store
deleted file mode 100644
Binary files a/src/.DS_Store and /dev/null differ
diff --git a/src/LazyPPL.hs b/src/LazyPPL.hs
--- a/src/LazyPPL.hs
+++ b/src/LazyPPL.hs
@@ -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'
 
diff --git a/src/LazyPPL/.DS_Store b/src/LazyPPL/.DS_Store
deleted file mode 100644
Binary files a/src/LazyPPL/.DS_Store and /dev/null differ
diff --git a/src/LazyPPL/Distributions.hs b/src/LazyPPL/Distributions.hs
--- a/src/LazyPPL/Distributions.hs
+++ b/src/LazyPPL/Distributions.hs
@@ -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)
diff --git a/src/LazyPPL/Distributions/.DS_Store b/src/LazyPPL/Distributions/.DS_Store
deleted file mode 100644
Binary files a/src/LazyPPL/Distributions/.DS_Store and /dev/null differ
diff --git a/src/LazyPPL/Distributions/IBP.hs b/src/LazyPPL/Distributions/IBP.hs
--- a/src/LazyPPL/Distributions/IBP.hs
+++ b/src/LazyPPL/Distributions/IBP.hs
@@ -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
diff --git a/src/LazyPPL/Distributions/Memoization.hs b/src/LazyPPL/Distributions/Memoization.hs
--- a/src/LazyPPL/Distributions/Memoization.hs
+++ b/src/LazyPPL/Distributions/Memoization.hs
@@ -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@ -}