diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -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 diff --git a/lazyppl.cabal b/lazyppl.cabal --- a/lazyppl.cabal +++ b/lazyppl.cabal @@ -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 . . @@ -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 diff --git a/src/LazyPPL.hs b/src/LazyPPL.hs --- a/src/LazyPPL.hs +++ b/src/LazyPPL.hs @@ -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. diff --git a/src/LazyPPL.hs-boot b/src/LazyPPL.hs-boot new file mode 100644 --- /dev/null +++ b/src/LazyPPL.hs-boot @@ -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 diff --git a/src/LazyPPL/Internal/ProbGen.hs b/src/LazyPPL/Internal/ProbGen.hs new file mode 100644 --- /dev/null +++ b/src/LazyPPL/Internal/ProbGen.hs @@ -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 #-} diff --git a/test/Api/Inference.hs b/test/Api/Inference.hs --- a/test/Api/Inference.hs +++ b/test/Api/Inference.hs @@ -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)