lazyppl 1.1.0 → 1.1.0.1
raw patch · 6 files changed
+172/−11 lines, 6 filesdep +mwc-randomPVP ok
version bump matches the API change (PVP)
Dependencies added: mwc-random
API changes (from Hackage documentation)
Files
- ChangeLog.md +7/−0
- lazyppl.cabal +13/−1
- src/LazyPPL.hs +10/−8
- src/LazyPPL.hs-boot +13/−0
- src/LazyPPL/Internal/ProbGen.hs +79/−0
- test/Api/Inference.hs +50/−2
ChangeLog.md view
@@ -1,5 +1,12 @@ # Changelog for lazyppl +## 1.1.0.1++### Fixed++- `wis`: when the log weights were large and equal, adding them rounded, so the cumulative weights all equalled the maximum and the first particle always won.+- `wis` now resamples in log space through the bounded categorical sampler of `mwc-random`, so equal large weights share the resamples; fixed-seed results change.+ ## 1.1.0 ### Added
lazyppl.cabal view
@@ -5,7 +5,7 @@ -- see: https://github.com/sol/hpack name: lazyppl-version: 1.1.0+version: 1.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/>. .@@ -50,6 +50,8 @@ LazyPPL.Internal.History LazyPPL.Internal.Sampler LazyPPL.Distributions.GP.Pure+ other-modules:+ LazyPPL.Internal.ProbGen hs-source-dirs: src build-depends:@@ -57,6 +59,7 @@ , containers >=0.6.7 && <0.9 , log-domain >=0.13.2 && <0.14 , math-functions >=0.3.4.2 && <0.4+ , mwc-random >=0.15.3.0 && <0.16 , random >=1.3.1 && <1.4 , transformers >=0.6.1.0 && <0.7 , vector >=0.13.2.0 && <0.14@@ -82,6 +85,7 @@ , lazyppl , log-domain >=0.13.2 && <0.14 , math-functions >=0.3.4.2 && <0.4+ , mwc-random >=0.15.3.0 && <0.16 , random >=1.3.1 && <1.4 , transformers >=0.6.1.0 && <0.7 , typed-process >=0.2.13.0 && <0.3@@ -103,6 +107,7 @@ , log-domain >=0.13.2 && <0.14 , math-functions >=0.3.4.2 && <0.4 , mtl >=2.3.1 && <2.4+ , mwc-random >=0.15.3.0 && <0.16 , random >=1.3.1 && <1.4 , transformers >=0.6.1.0 && <0.7 , vector >=0.13.2.0 && <0.14@@ -120,6 +125,7 @@ , lazyppl , log-domain >=0.13.2 && <0.14 , math-functions >=0.3.4.2 && <0.4+ , mwc-random >=0.15.3.0 && <0.16 , random >=1.3.1 && <1.4 , transformers >=0.6.1.0 && <0.7 , vector >=0.13.2.0 && <0.14@@ -137,6 +143,7 @@ , lazyppl , log-domain >=0.13.2 && <0.14 , math-functions >=0.3.4.2 && <0.4+ , mwc-random >=0.15.3.0 && <0.16 , random >=1.3.1 && <1.4 , transformers >=0.6.1.0 && <0.7 , vector >=0.13.2.0 && <0.14@@ -156,6 +163,7 @@ , lazyppl , log-domain >=0.13.2 && <0.14 , math-functions >=0.3.4.2 && <0.4+ , mwc-random >=0.15.3.0 && <0.16 , random >=1.3.1 && <1.4 , transformers >=0.6.1.0 && <0.7 , vector >=0.13.2.0 && <0.14@@ -175,6 +183,7 @@ , lazyppl , log-domain >=0.13.2 && <0.14 , math-functions >=0.3.4.2 && <0.4+ , mwc-random >=0.15.3.0 && <0.16 , random >=1.3.1 && <1.4 , transformers >=0.6.1.0 && <0.7 , vector >=0.13.2.0 && <0.14@@ -191,6 +200,7 @@ , lazyppl , log-domain >=0.13.2 && <0.14 , math-functions >=0.3.4.2 && <0.4+ , mwc-random >=0.15.3.0 && <0.16 , random >=1.3.1 && <1.4 , transformers >=0.6.1.0 && <0.7 , vector >=0.13.2.0 && <0.14@@ -209,6 +219,7 @@ , lazyppl , log-domain >=0.13.2 && <0.14 , math-functions >=0.3.4.2 && <0.4+ , mwc-random >=0.15.3.0 && <0.16 , random >=1.3.1 && <1.4 , transformers >=0.6.1.0 && <0.7 , vector >=0.13.2.0 && <0.14@@ -226,6 +237,7 @@ , lazyppl , log-domain >=0.13.2 && <0.14 , math-functions >=0.3.4.2 && <0.4+ , mwc-random >=0.15.3.0 && <0.16 , random >=1.3.1 && <1.4 , transformers >=0.6.1.0 && <0.7 , vector >=0.13.2.0 && <0.14
src/LazyPPL.hs view
@@ -61,10 +61,14 @@ import Control.Monad.Trans.Writer import Control.Monad.Trans.Class import Data.Monoid+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU import System.Random hiding (uniform)+import qualified System.Random.MWC.Distributions as MWC import Control.Monad import System.IO.Unsafe (unsafeInterleaveIO) import LazyPPL.Internal.History (History, newHistory, appendMutation, labelAt)+import LazyPPL.Internal.ProbGen (ProbGen(..)) import LazyPPL.Internal.Sampler (randomLabel, checkedLogWeight) import Numeric.Log @@ -179,21 +183,19 @@ {- | Weighted importance sampling first draws n weighted samples, and then samples a stream of results from that, regarded as an empirical distribution. Sometimes called "likelihood weighted importance sampling". -This is a reference implementation. It will not usually be very efficient at all, but may be useful for debugging. Raises an error if all @n@ samples have weight 0.+This is a reference implementation. It will not usually be very efficient at all, but may be useful for debugging. Raises an error if @n@ is less than 1 or all @n@ samples have weight 0. -} wis :: Int -- ^ @n@, the number of samples to base on, at least 1 -> Meas a -- ^ @m@, the measure to normalize -> IO [a] -- ^ Returns a stream of samples+wis n _ | n <= 0 = error "wis: the number of samples must be positive" wis n m = do xws <- weightedsamples m- let xws' = take n $ accumulate xws 0- let max = snd $ last xws'+ let particles = V.fromList (take n xws)+ weights = VU.fromList (map (ln . snd) (V.toList particles)) g <- newStdGen- let rs = (randoms g :: [Double])- if max == 0 then error "wis: every sample has weight 0"- else return $ map (\r -> fst $ head $ filter (\(x, w) -> w > 0 && w >= Exp (log r) * max) xws') rs- where accumulate ((x, w) : xws) a = (x, w + a) : accumulate xws (w + a)- accumulate [] a = []+ if VU.all (== (-1 / 0)) weights then error "wis: every sample has weight 0"+ else return $ map (\tree -> fst (particles V.! runProb (MWC.logCategorical weights ProbGen) tree)) (randomTrees g) {- | Produce a stream of samples, using Metropolis Hastings simulation.
+ src/LazyPPL.hs-boot view
@@ -0,0 +1,13 @@+module LazyPPL where++import Control.Monad.Trans.Writer+import Data.Monoid+import Numeric.Log++data Tree = Tree Double [Tree]+newtype Prob a = Prob (Tree -> a)+newtype Meas a = Meas (WriterT (Product (Log Double)) Prob a)++instance Monad Prob+instance Functor Prob+instance Applicative Prob
+ src/LazyPPL/Internal/ProbGen.hs view
@@ -0,0 +1,79 @@+{-# LANGUAGE MultiParamTypeClasses #-}++{- | The stock 'StatefulGen' view of 'LazyPPL.Prob', used by 'LazyPPL.wis' to+resample through the bounded samplers of @mwc-random@. This module is internal+and not exposed, so the 1.1.0 interface of the exposed modules is unchanged:+the instance below mentions only the internal 'ProbGen' and the exported+'Prob', and haddock lists it under this hidden module, not under 'LazyPPL'. -}+module LazyPPL.Internal.ProbGen (ProbGen(..)) where++import Control.Monad.ST (runST)+import Data.Bits ((.|.), shiftL, shiftR)+import Data.STRef (newSTRef, readSTRef, writeSTRef)+import Data.Word (Word64)+import {-# SOURCE #-} LazyPPL (Prob(..), Tree(..))+import System.Random.Stateful (StatefulGen(..), fillByteArrayST)++{- | A stateful generator view of a probability computation, for the stock+sampling algorithms of @random@ and @mwc-random@ that take a 'StatefulGen',+such as 'System.Random.MWC.Distributions.logCategorical'.++Precondition: the tree must provide a child for every requested word; the+library's trees always do (a hand-built tree that runs out of children raises+\"LazyPPL: random tree exhausted\").++'LazyPPL.Internal.Sampler.randomLabel' puts every label on the midpoint grid+(k + 1/2) / 2^52 with k a uniform 52-bit integer, and 'labelIndex' recovers k+exactly, so each label supplies 52 uniform bits. 'uniformWord64' joins the+label at its tree node with the label of the node's first child — a subtree+that 'Prob' sequencing reserves to this primitive — into 64 uniform bits; the+smaller word methods use one label. Every request a stock algorithm makes,+including the repeated draws of its rejection loops, reads a distinct tree+coordinate under the monad structure of 'Prob', and running the same+computation on the same tree gives the same words. -}+data ProbGen = ProbGen++-- The 52-bit grid index of a label, the k of the midpoint grid+-- u = (k + 1/2) / 2^52 of 'LazyPPL.Internal.Sampler.randomLabel'. It is exact+-- for every label 'randomLabel' draws: u * 2^52 = k + 1/2 is exactly+-- representable in binary64, and floor (k + 1/2) is k. Labels of hand-built+-- trees keep the index of the grid point at or below them.+labelIndex :: Double -> Word64+labelIndex u = fromInteger (floor (u * 4503599627370496))++instance StatefulGen ProbGen Prob where+ uniformWord64 ProbGen = Prob $ \(Tree r ts) -> case ts of+ (Tree c _ : _) -> (labelIndex r `shiftL` 12) .|. (labelIndex c `shiftR` 40)+ [] -> treeExhausted+ {-# INLINE uniformWord64 #-}+ uniformWord32 ProbGen = Prob $ \(Tree r _) -> fromIntegral (labelIndex r `shiftR` 20)+ {-# INLINE uniformWord32 #-}+ -- The class default delegates to a frozen generator (System.Random.Internal+ -- 1.3.1, uniformByteArrayM's default signature), which a tree coordinate has+ -- no notion of. The stock fillByteArrayST (System.Random.Stateful, under+ -- "Helper functions for creating instances") fills the array little-endian+ -- from successive words, drawn here from successive children through an+ -- STRef; the empty list is matched explicitly, so no byte stays+ -- uninitialised.+ uniformByteArrayM isPinned n ProbGen = Prob $ \(Tree _ gs) -> runST $ do+ ref <- newSTRef gs+ fillByteArrayST isPinned n $ do+ ts <- readSTRef ref+ case ts of+ (t : ts') -> do+ writeSTRef ref ts'+ pure (treeWord t)+ [] -> treeExhausted+ {-# INLINE uniformByteArrayM #-}++-- One clear error when a hand-built tree runs out of children.+treeExhausted :: a+treeExhausted = error "LazyPPL: random tree exhausted"+{-# NOINLINE treeExhausted #-}++-- One tree word per node, as uniformWord64 makes it.+treeWord :: Tree -> Word64+treeWord (Tree r ts) = case ts of+ (Tree c _ : _) -> (labelIndex r `shiftL` 12) .|. (labelIndex c `shiftR` 40)+ [] -> treeExhausted+{-# INLINE treeWord #-}
test/Api/Inference.hs view
@@ -64,7 +64,30 @@ pure x) positiveSupport <- assertBool "wis returns supported samples from a batch containing zero weights" (all (>0.5) positive)- pure (escaped : retained : moving ++ [noZeroNormalization, positiveSupport])+ emptyBatch <- try (wis 0 (pure ()) >>= evaluate . take 1) :: IO (Either ErrorCall [()])+ emptyRejected <- assertBool "wis rejects an empty empirical measure"+ (either (const True) (const False) emptyBatch)+ let base = do x <- sample uniform; scoreLog (Exp 0); pure x+ setStdGen (mkStdGen 29)+ [(target, _)] <- take 1 <$> weightedsamples base+ let oneModel = do+ x <- sample uniform+ scoreLog (Exp (if x == target then 0 else -1 / 0))+ pure x+ setStdGen (mkStdGen 29)+ oneParticles <- take 64 <$> weightedsamples oneModel+ setStdGen (mkStdGen 29)+ oneResamples <- take 2000 <$> wis 64 oneModel+ onePositive <- assertBool "the batch has one positive particle and only it is resampled"+ (length (filter ((> 0) . snd) oneParticles) == 1 && all (== target) oneResamples)+ let sepModel = do+ x <- sample (bernoulli 0.5)+ scoreLog (Exp (if x then -1000 else 0))+ pure x+ setStdGen (mkStdGen 24)+ sepResamples <- take 2000 <$> wis 64 sepModel+ separated <- assertBool "negligible-weight particles are not sampled" (all not sepResamples)+ pure (escaped : retained : moving ++ [noZeroNormalization, positiveSupport, emptyRejected, onePositive, separated]) -- A fair coin b. On true, a uniform u is read and scored 6u, which integrates -- to 3; on false, the score is 1. So P(b) = 3/4 exactly, and the site u exists@@ -110,10 +133,35 @@ assertBool ("wis 5 resamples each particle in proportion to its weight, (observed, expected) = " ++ show shares) (all within5 shares && length (filter (`elem` map fst particles) resamples) == draws) +-- Resampling in log space: equal weights share the resamples equally no+-- matter how large the log weights are. With log weight 1e308 the cumulative+-- sums of the weights all round to the maximum, so a resampler that adds the+-- weights always picks the first particle; this regression check observes the+-- share of one coin instead.+wisLogSpaceChecks :: IO [Bool]+wisLogSpaceChecks = do+ let equalModel = do+ x <- sample (bernoulli 0.5)+ scoreLog (Exp 1e308)+ pure x+ n = 64+ draws = 4000+ setStdGen (mkStdGen 29)+ equalParticles <- take n <$> weightedsamples equalModel+ setStdGen (mkStdGen 29)+ equalResamples <- take draws <$> wis n equalModel+ let expected = fromIntegral (length (filter fst equalParticles)) / fromIntegral n :: Double+ observed = fromIntegral (length (filter id equalResamples)) / fromIntegral draws+ equalLarge <- assertBool ("equal-weight wis share, expected " ++ show expected ++ ", got " ++ show observed)+ (expected > 0 && expected < 1 && abs (observed - expected) < 0.05)+ pure [equalLarge]+ posteriorChecks :: IO [Bool] posteriorChecks = do setStdGen (mkStdGen 21) chain <- map fst . take 40000 . drop 1000 <$> mh 0.5 coinModel setStdGen (mkStdGen 22) irreducible <- map fst . take 40000 . drop 1000 <$> mhirreducible 0.5 0.1 coinModel- sequence [posteriorFraction "mh 0.5 coin posterior" chain, posteriorFraction "mhirreducible 0.5 0.1 coin posterior" irreducible, wisResamplingCheck]+ logSpace <- wisLogSpaceChecks+ rest <- sequence [posteriorFraction "mh 0.5 coin posterior" chain, posteriorFraction "mhirreducible 0.5 0.1 coin posterior" irreducible, wisResamplingCheck]+ pure (logSpace ++ rest)