mcmc-samplers 0.1.0.0 → 0.1.1.0
raw patch · 14 files changed
+1708/−528 lines, 14 filesdep +hakarudep −vectordep ~basedep ~hmatrixdep ~statistics
Dependencies added: hakaru
Dependencies removed: vector
Dependency ranges changed: base, hmatrix, statistics
Files
- Actions.hs +0/−72
- Distributions.hs +0/−257
- Kernels.hs +0/−103
- MCMC/Actions.hs +65/−0
- MCMC/Combinators.hs +206/−0
- MCMC/Distributions.hs +215/−0
- MCMC/Examples/GMM.hs +475/−0
- MCMC/Examples/HandwrittenGMM.hs +311/−0
- MCMC/Kernels.hs +111/−0
- MCMC/SemanticEditors.hs +74/−0
- MCMC/Tests.hs +125/−0
- MCMC/Types.hs +110/−0
- Tests.hs +0/−88
- mcmc-samplers.cabal +16/−8
− Actions.hs
@@ -1,72 +0,0 @@-module Actions ( Action (..)- , execute- , every- , Batch- , BatchAct- , BatchAction- , inBatches- , pack- , PrintF- , batchViz- , batchPrint- ) where--import Control.Monad--type Act x m a = x -> a -> m a--data Action x m a b = Action (Act x m a) (a -> m b) a--execute :: Monad m => Action x m a b -> x -> m (Action x m a b)-execute (Action act fin a) x = liftM (Action act fin) (act x a)--every :: Monad m => Int -> Action x m a b -> Action x m (a,Int) b-every n (Action act fin a) = - let skip_act x (b,i) = if i == (n-1)- then do b' <- act x b- return (b',0)- else return (b,i+1)- skip_fin = fin . fst- in Action skip_act skip_fin (a,0)---- Batch actions -- --type Batch x = ([x], Int)-type BatchAct x m = Act x m (Batch x)-type BatchAction x m b = Action x m (Batch x) b--inBatches :: Monad m => (Batch x -> m b) -> Int -> BatchAct x m-inBatches f n x a@(l, i)- | i == n = f a >> return ([x], 1)- | otherwise = return (x:l, i+1)--pack :: (Batch x -> m b) -> BatchAct x m -> BatchAction x m b-pack f act = Action act f ([], 0)---- Visualization ----type PrintF x s = [x] -> [s]--vizJSON :: Show s => [s] -> String-vizJSON samplelist = "{\"rvars\": {\"x\": " ++ show samplelist ++ "}}"--closeJSON :: Show s => [s] -> String-closeJSON samplelist = "{\"close\": true, \"rvars\": {\"x\": " ++ show samplelist ++ "}}"--visualize :: Show s => PrintF x s -> Batch x -> IO ()-visualize f (ls,_) = unless (null ls) $ putStrLn.vizJSON $ f ls--vizClose :: Show s => PrintF x s -> Batch x -> IO ()-vizClose f (ls,_) = unless (null ls) $ putStrLn.closeJSON $ f ls--batchViz :: Show s => PrintF x s -> Int -> BatchAction x IO ()-batchViz f n = - let viz = visualize f - close = vizClose f- in pack close $ inBatches viz n--batchPrint :: Show s => PrintF x s -> Int -> BatchAction x IO ()-batchPrint f n =- let p fn (ls,_) = unless (null ls) $ print $ fn ls- in pack (p f) $ inBatches (p f) n-
− Distributions.hs
@@ -1,257 +0,0 @@-{-# LANGUAGE MultiParamTypeClasses, KindSignatures, FlexibleInstances, GADTs, - OverlappingInstances #-}--module Distributions ( Rand- , Probability- , Target (..)- , makeTarget- , density- , Proposal (..)- , makeProposal- , sampleFrom- , fromProposal- , uniform- , diag- , normal- , categorical- , targetMix- , proposalMix- , first- , second- , car- , cdr- , nth- , block- , swapWith- , updateNth- , updateBlock- ) where--import qualified System.Random.MWC as MWC-import qualified System.Random.MWC.Distributions as MWC.D-import Control.Monad-import Control.Monad.Primitive-import qualified Data.Packed.Matrix as M-import qualified Numeric.LinearAlgebra.Algorithms as LA-import qualified Numeric.Container as C-import qualified Data.Vector as V--type Rand = MWC.Gen (PrimState IO)-type Probability = Double--type Density a = a -> Probability--class HasDensity d a where- density :: d a -> Density a--data Target a = T (Density a)--makeTarget :: Density a -> Target a-makeTarget = T--instance HasDensity Target a where- density (T d) = d--type Sample a = Rand -> IO a-data Proposal a = P (Density a) (Sample a)--makeProposal :: Density a -> Sample a -> Proposal a-makeProposal = P--instance HasDensity Proposal a where- density (P d _) = d--sampleFrom :: Proposal a -> Sample a-sampleFrom (P _ s) = s--fromProposal :: Proposal a -> Target a-fromProposal = T . density---- Uniform -- --uniform :: (MWC.Variate a, Real a) => [a] -> [a] -> Proposal [a]-uniform a b- | b < a = uniform b a- | a < b = makeUniform a b- | otherwise = error "Wrong parameters for Uniform distribution"--unif1D :: Real a => a -> a -> a -> Probability-unif1D a b x- | x < a = 0- | x > b = 0- | otherwise = 1 / realToFrac (b - a)--makeUniform :: (MWC.Variate a, Real a) => [a] -> [a] -> Proposal [a]-makeUniform a b = - let tuf f (p,q,r) = f p q r- uniD x = product . map (tuf unif1D) $ zip3 a b x- uniSF g = mapM (flip MWC.uniformR g) $ zip a b- in P uniD uniSF- --- Normal ----type CovMatrix = M.Matrix Double-type Mu a = M.Matrix a--mu :: M.Element a => [a] -> Mu a-mu mean = M.fromLists [mean]--diag :: [Double] -> [[Double]]-diag d = [(nth i) (swapWith e) (replicate (length d) 0) | (i,e) <- zip [1..] d]--normal :: [Double] -> [[Double]] -> Proposal [Double]-normal mean cov =- let covMat = M.fromLists cov- (muMat, n) = (mu mean, length mean)- in P (normalD muMat covMat) (normalSF muMat covMat n)--normalD :: Mu Double -> CovMatrix -> Density [Double]-normalD m cov x = c * exp (-d / 2)- where (covInv, (lndet, sign)) = LA.invlndet cov- c1 = (2*pi) ^^ (length x)- c = 1 / (sqrt $ sign * (exp lndet) * c1)- xm = C.sub (M.fromLists [x]) m- prod = xm C.<> covInv C.<> (M.trans xm)- d = (M.@@>) prod (0,0)--normalSF :: Mu Double -> CovMatrix -> Int -> Sample [Double]-normalSF m cov n g = do- z <- replicateM n (MWC.D.standard g)- let zt = M.trans $ M.fromLists [z]- a = LA.chol cov- return . head . M.toLists $ C.add m $ C.trans $ a C.<> zt---- Categorical ----categorical :: [a] -> [Probability] -> Proposal (a,Int)-categorical cs ps = let (cats,probs) = vectorCats cs ps in makeCategorical cats probs---- The index is stored (a,Int) to allow categorical dists over objects outside the Eq class--- Assumption: --- 1. length of probability list >= length of category list - 1--- 2. the order in each list "matters", i.e, first prob maps to first category, second to second, and so on-makeCategorical :: V.Vector a -> V.Vector Probability -> Proposal (a,Int)-makeCategorical cats probs = let den (_,i) = probs V.! i- in P den (catSF cats probs)--catSF :: V.Vector a -> V.Vector Probability -> Sample (a,Int)-catSF cats probs g = do - u <- sampleFrom (uniform [0] [1]) g- let cdfs = V.prescanl (+) 0 probs- i = V.maxIndex $ V.filter ((>=) $ head u) cdfs- return $ (cats V.! i, i)--normalize :: V.Vector Probability -> V.Vector Probability-normalize probs = let s = V.sum probs in V.map (flip (/) s) probs--vectorCats :: [a] -> [Probability] -> (V.Vector a, V.Vector Probability)-vectorCats cs ps = (V.fromList cs, normalize $ V.fromList ps)---- Target Mixtures ----mixD :: HasDensity t a => V.Vector (t a) -> V.Vector Probability -> Density a-mixD cats probs x = V.foldl1 (+) $ V.imap (\i t -> (probs V.! i)*(density t x)) cats--targetMix :: HasDensity t a => [t a] -> [Probability] -> Target a-targetMix ts ps =- let (cats,probs) = vectorCats ts ps- in T (mixD cats probs)- --- Proposal Mixtures ----proposalMix :: [Proposal a] -> [Probability] -> Proposal a-proposalMix props ps = - let (cats,probs) = vectorCats props ps- propCat = makeCategorical cats probs- mixSF g = do- (prop,_) <- sampleFrom propCat g- sampleFrom prop g- in P (mixD cats probs) mixSF---- Semantic editor combinators ------ http://conal.net/blog/posts/semantic-editor-combinators--first :: (a -> a') -> ((a,b) -> (a',b))-second :: (b -> b') -> ((a,b) -> (a,b'))- -first f = \ (a,b) -> (f a, b)-second g = \ (a,b) -> (a, g b)--car :: (a -> a) -> ([a] -> [a])-cdr :: ([a] -> [a]) -> ([a] -> [a])--car f = \(x:xs) -> f x : xs-cdr f = \(x:xs) -> x : f xs--nth :: Int -> (a -> a) -> ([a] -> [a])-nth 1 = car-nth n = cdr . nth (n-1)--carM :: Monad m => (a -> m a) -> ([a] -> m [a])-carM f (x:xs) = do x' <- f x - return $ x' : xs--cdrM :: Monad m => ([a] -> m [a]) -> ([a] -> m [a])-cdrM f (x:xs) = do xs' <- f xs- return $ x : xs'--nthM :: Monad m => Int -> (a -> m a) -> ([a] -> m [a])-nthM 1 = carM-nthM n = cdrM . nthM (n-1)--block :: Int -> Int -> ([a] -> [a]) -> ([a] -> [a])-block begin end f ls- | begin == end = nth begin (unlift f) ls- | begin == 1 = f (take end ls) ++ drop end ls- | otherwise = front ++ f mids ++ back- where (front, mids, back) = chopAt begin end ls--blockM :: Monad m => Int -> Int -> ([a] -> m [a]) -> ([a] -> m [a])-blockM begin end f ls- | begin == end = nthM begin (unliftM f) ls- | begin == 1 = f (take end ls) >>= return . flip (++) (drop end ls)- | otherwise = do let (front, mids, rest) = chopAt begin end ls- fmids <- f mids- return $ front ++ fmids ++ rest--swapWith :: a -> (b -> a)-swapWith x _ = x--unlift :: ([a] -> [a]) -> a -> a-unlift f x = head $ f [x]--unliftM :: Monad m => ([a] -> m [a]) -> a -> m a-unliftM f x = f [x] >>= return.head---- A sample is 1-indexed, i.e., dimensions go from 1 to n-getBlock :: Int -> Int -> [a] -> [a]-getBlock begin end- | begin == end = \ls -> [ls !! (begin - 1)]- | otherwise = take (end + 1 - begin) . drop (begin-1)--chopAt :: Int -> Int -> [a] -> ([a], [a], [a])-chopAt begin end ls = (front, mids, back)- where (front, rest) = splitAt (begin-1) ls- (mids, back) = splitAt (end + 1 - begin) rest--updateNth :: Int -> ([a] -> Proposal [a]) -> ([a] -> Proposal [a])-updateNth n p x = - let den y = flip density (getBlock n n y) $ p (getBlock n n x)- s g = nthM n (unliftM (\xn -> sampleFrom (p xn) g)) x- in P den s--updateBlock :: Int -> Int -> ([a] -> Proposal [a]) -> [a] -> Proposal [a]-updateBlock n m p x = - let den y = flip density (getBlock n m y) $ p (getBlock n m x)- s g = blockM n m (\b -> sampleFrom (p b) g) x- in P den s---- ex = (second.first.second) not (1,((3,True),2))--- eg = (cdr.cdr.car) not [False,False,True,False]--- examp = (nth 4) not [False,False,False,True,False]--- exampl = (nth 4) (swapWith 4) [1,2,3,5,5]--- diagex = diag [1,2,3,4]--- cdrex = cdr (\ls -> map ((+)10) ls) [1,2,3,4,5]--- cprobs = U.prescanl (+) 0 $ normalize $ U.fromList [0.3, 0.4, 0.6, 0.3, 0.7] --- gcp = U.maximum $ U.filter ((>=) 0.60) cprobs
− Kernels.hs
@@ -1,103 +0,0 @@-{-# LANGUAGE GADTs, MultiParamTypeClasses, KindSignatures, FlexibleInstances #-}--module Kernels ( Step- , Kernel- , walk- , metropolisHastings- , vizMH- , printMH- , Temp- , CoolingSchedule- , StateSA- , simulatedAnnealing- , vizSA- , printSA- , mixSteps- , cycleKernel- ) where--import Distributions-import Actions-import Text.Printf---- Kernels ----type Step x = Rand -> x -> IO x-type Kernel x a = Target a -> (a -> Proposal a) -> Step x-walk :: Step x -> x -> Int -> Rand -> Action x IO a b -> IO b-walk _ _ 0 _ (Action _ f a) = f a-walk step x n r action = do - x' <- step r x- execute action x' >>= walk step x' (n-1) r---- Metropolis Hastings ----metropolisHastings :: Kernel [a] [a]-metropolisHastings t c_p = - let mhStep g xi = do- u <- sampleFrom (uniform [0] [1]) g- xstar <- sampleFrom (c_p xi) g- let accept = min 1 (numer / denom)- numer = density t xstar * density (c_p xstar) xi- denom = density t xi * density (c_p xi) xstar- return $ if head u < accept then xstar else xi- in mhStep---- Visualizes only the first dimension-vizMH :: PrintF [Double] Double-vizMH = map head--printMH :: PrintF [Double] [String]-printMH lls = let p d = printf "%0.3f" d :: String- in map (map p) lls---- Simulated Annealing ----type Temp = Double-type CoolingSchedule = Temp -> Temp-type StateSA a = (a, Temp, CoolingSchedule)--simulatedAnnealing :: Kernel (StateSA [a]) [a]-simulatedAnnealing t c_p = - let saStep g (xi,temp,cool) = do- u <- sampleFrom (uniform [0] [1]) g- xstar <- sampleFrom (c_p xi) g- let accept = min 1 (numer / denom)- numer = (*) (density (c_p xstar) xi) $ (**) (1 / temp) (density t xstar)- denom = (*) (density (c_p xi) xstar) $ (**) (1 / temp) (density t xi)- new_temp = cool temp- return $ if head u < accept then (xstar,new_temp,cool) else (xi,new_temp,cool)- in saStep--tripleFirst :: (a, b, c) -> a-tripleFirst (a,_,_) = a--myFilter :: [[Double]] -> [[Double]]-myFilter = filter (\x -> x < (repeat 15) && x > (repeat $ -5))---- Visualizes only the first dimension-vizSA :: PrintF (StateSA [Double]) Double-vizSA = vizMH . myFilter . map tripleFirst---- Print the samples without the temp/cooling schedule-printSA :: PrintF (StateSA [Double]) [String]-printSA = printMH . myFilter . map tripleFirst---- Kernel Mixtures ----mixSteps :: [Step x] -> [Probability] -> Step x-mixSteps steps probs = - let mixStep g x = do- (step,_) <- sampleFrom (categorical steps probs) g- step g x- in mixStep---- Kernel Cycles ----cycleKernel :: Kernel x a -> Target a -> [a -> Proposal a] -> Step x-cycleKernel kernel t cps =- let steps g = [kernel t cp g | cp <- cps]- combine comb step = (\iox -> iox >>= comb) . step- cycleStep g = foldl combine return (steps g)- in cycleStep-
+ MCMC/Actions.hs view
@@ -0,0 +1,65 @@+module MCMC.Actions ( Action (..)+ -- * Predefined actions+ , collect+ , display+ -- ** Batch actions+ , Batch+ , BatchAct+ , BatchAction+ , inBatches+ , pack+ , PrintF+ , batchViz+ , batchPrint+ ) where++import Control.Monad+import MCMC.Types++collect :: Action x IO [x] [x]+collect = makeAction (\ x ls -> return (x:ls)) return []++display :: Show s => (x -> s) -> Action x IO () ()+display f = makeAction (\x _ -> print $ f x) print ()++-- Batch actions -- ++type Batch x = ([x], Int)+type BatchAct x m = Act x m (Batch x)+type BatchAction x m b = Action x m (Batch x) b++inBatches :: Monad m => (Batch x -> m b) -> Int -> BatchAct x m+inBatches f n x a@(l, i)+ | i == n = f a >> return ([x], 1)+ | otherwise = return (x:l, i+1)++pack :: (Batch x -> m b) -> BatchAct x m -> BatchAction x m b+pack f act = makeAction act f ([], 0)++-- Batch Visualization --++type PrintF x s = [x] -> [s]++vizJSON :: Show s => [s] -> String+vizJSON samplelist = "{\"rvars\": {\"x\": " ++ show samplelist ++ "}}"++closeJSON :: Show s => [s] -> String+closeJSON samplelist = "{\"close\": true, \"rvars\": {\"x\": " ++ show samplelist ++ "}}"++visualize :: Show s => PrintF x s -> Batch x -> IO ()+visualize f (ls,_) = unless (null ls) $ putStrLn.vizJSON $ f ls++vizClose :: Show s => PrintF x s -> Batch x -> IO ()+vizClose f (ls,_) = unless (null ls) $ putStrLn.closeJSON $ f ls++batchViz :: Show s => PrintF x s -> Int -> BatchAction x IO ()+batchViz f n = + let viz = visualize f + close = vizClose f+ in pack close $ inBatches viz n++batchPrint :: Show s => PrintF x s -> Int -> BatchAction x IO ()+batchPrint f n =+ let p fn (ls,_) = unless (null ls) $ print $ fn ls+ in pack (p f) $ inBatches (p f) n+
+ MCMC/Combinators.hs view
@@ -0,0 +1,206 @@+{-# LANGUAGE ViewPatterns #-}++module MCMC.Combinators (+ -- * Target combinators+ mixTargets+ -- * Proposal combinators+ , mixProposals+ , chooseProposal+ , mixCondProposals+ , updateNth+ , updateBlock+ , updateFirst+ , updateSecond+ -- * Step combinators+ , mixSteps+ , chooseStep+ , cycleStep+ -- * Action combinators+ , execute+ , every+ ) where++import MCMC.Types+import MCMC.Distributions+import MCMC.SemanticEditors+import Data.Maybe++import Control.Monad++-- Target Mixtures --++-- Assume normalized probabilities+mixDensity :: HasDensity t a => [(t a, Double)] -> Density a+mixDensity catProbs x = foldl1 (+) $ map f catProbs+ where f (t,p) = p*(density t x)++-- | Create a mixture target distribution from a list of +-- distribution-weight pairs. +-- +-- The input distributions can be of any type in the 'HasDensity'+-- class, which includes 'Proposal' and 'Target'.+-- The output is a 'Target' containing a 'Density'.+-- +-- The input weights are meant to be relative, not+-- required to be normalized.+mixTargets :: HasDensity t a => [(t a, Double)] -> Target a+mixTargets = makeTarget . mixDensity . normalize+ +-- Proposal Mixtures --++-- | Create a mixture proposal distribution from a list of+-- distribution-weight pairs.+-- +-- The input weights are meant to be relative, not+-- required to be normalized.+mixProposals :: [(Proposal a, Double)] -> Proposal a+mixProposals proposalProbs = + let normedPropProbs = normalize proposalProbs+ mixSF g = do+ let (props, probs) = unzip normedPropProbs+ propMap = zip [1..] props+ propKey <- sampleFrom (categoricalNormed (zip [1..] probs)) g+ let proposal = fromMaybe (error "mixProposals: undefined key")+ (lookup propKey propMap)+ sampleFrom proposal g+ in makeProposal (mixDensity normedPropProbs) mixSF++-- Proposal Choices -- ++-- | Create a uniform mixture of proposals based on the mapping+-- function. The first argument is the total number of proposals+-- in the mixture.+-- +-- This is a convenience function for cases where there are a large +-- number of index-based proposals (such as the case where there is a +-- unique proposal for updating each dimension in a high-dimensional state).+chooseProposal :: Int -> (Int -> Proposal a) -> Proposal a+chooseProposal n f = mixProposals $ [(f i , 1) | i <- [1..n]] -- uniform mix++-- Proposal updaters combinations++-- | Create a mixture of /conditional/ proposal distributions.+-- The result is a mixture of proposal distributions that are all conditioned+-- on the same starting state.+-- +-- The input weights are meant to be relative, not+-- required to be normalized.+mixCondProposals :: [(a -> Proposal a, Double)] -> a -> Proposal a+mixCondProposals cps a = mixProposals $ map (\(cp,prob) -> (cp a, prob)) cps++getNth :: Int -> [a] -> ([a], a, [a])+getNth n ls = let (front, [e], back) = chopAt n n ls+ in (front, e, back)++-- | Update the nth dimension of a multi-dimensioanl state +-- represented using a list.+updateNth :: Eq a => Int -- ^ The dimension index+ -> (a -> Proposal a) -- ^ A conditional proposal to focus on that dimension+ -> ([a] -> Proposal [a]) -- ^ A conditional proposal that updates only one dimension+updateNth n p x = + let den y = let (xBefore, xn, xAfter) = getNth n x+ (yBefore, yn, yAfter) = getNth n y+ in if length x == length y && xBefore == yBefore+ && xAfter == yAfter+ then density (p xn) yn else 0+ s g = nthM n (\xn -> sampleFrom (p xn) g) x+ in makeProposal den s++-- | Update a contiguous block of dimensions of a multi-dimensional state+-- represented using a list.+updateBlock :: Eq a => Int -- ^ The start dimension of the block + -> Int -- ^ The end dimension of the block+ -> ([a] -> Proposal [a]) -- ^ A conditional proposal to focus on the block + -> ([a] -> Proposal [a]) -- ^ A conditional proposal that updates only that block+updateBlock n m p x = + let den y = let (xBefore, xBlock, xAfter) = chopAt n m x+ (yBefore, yBlock, yAfter) = chopAt n m y+ in if length x == length y && xBefore == yBefore+ && xAfter == yAfter+ then density (p xBlock) yBlock else 0+ s g = blockM n m (\b -> sampleFrom (p b) g) x+ in makeProposal den s++-- | Update the first element of a multi-dimensional state represented as a tuple.+updateFirst :: Eq b => (a -> Proposal a) -- ^ A conditional proposal to focus on the first element + -> ((a,b) -> Proposal (a,b))+updateFirst p (a,b) = makeProposal dens sf+ where dens (x,y) = if y==b then density (p a) x else 0+ sf g = do a' <- sampleFrom (p a) g+ return (a',b)++-- | Update the second element of a multi-dimensional state represented as a tuple.+updateSecond :: Eq a => (b -> Proposal b) + -> ((a,b) -> Proposal (a,b)) -- ^ A conditional proposal to focus on the second element+updateSecond p (a,b) = makeProposal dens sf+ where dens (x,y) = if x==a then density (p b) y else 0+ sf g = do b' <- sampleFrom (p b) g+ return (a,b')++-- Examples --++-- ex = (second.first.second) not (1,((3,True),2))+-- eg = (cdr.cdr.car) not [False,False,True,False]+-- examp = (nth 4) not [False,False,False,True,False]+-- exampl = (nth 4) (swapWith 4) [1,2,3,5,5]+-- diagex = diag [1,2,3,4]+-- cdrex = cdr (\ls -> map ((+)10) ls) [1,2,3,4,5]+-- cprobs = U.prescanl (+) 0 $ normalize $ U.fromList [0.3, 0.4, 0.6, 0.3, 0.7]+-- gcp = U.maximum $ U.filter ((>=) 0.60) cprobs++-- Kernel Mixtures --++-- | Create a mixture 'Step', representing a categorical distribution over +-- multiple inference procedures.+-- +-- The input weights are meant to be relative, not+-- required to be normalized.+mixSteps :: [(Step x, Double)] -> Step x+mixSteps stepProbs = + let mixStep g x = do+ let (steps, probs) = unzip stepProbs+ stepMap = zip [1..] steps+ stepKey <- sampleFrom (categorical (zip [1..] probs)) g+ let step = fromMaybe (error "mixSteps: undefined key")+ (lookup stepKey stepMap)+ step g x+ in mixStep++-- | The analogue of 'chooseProposal', to use a specific 'Step'+-- for each dimension of a multi-dimensional state.+-- +-- The same 'Kernel' value is used for creating each 'Step' in the mixture.+chooseStep :: Int -- ^ The total number of 'Step's in the mixture+ -> (Int -> Target a) -- ^ The mapping function over targets+ -> (Int -> a -> Proposal a) -- ^ The mapping function over conditional proposals+ -> Kernel x a -- ^ The transition kernel to use across all 'Step's+ -> Step x -- ^ The mixture step+chooseStep n t p k = mixSteps [(k (t i) (p i) , 1) | i <- [1..n]]++-- Kernel Cycles --++-- | Create a cycled transition kernel. This combinator applies the +-- provided transition kernel to the target, and to the conditional proposals+-- in the order that they appear in the list.+-- Each application generates an intermediate state that is passed to the next +-- conditional proposal in the list.+cycleStep :: Kernel x a -> Target a -> [a -> Proposal a] -> Step x+cycleStep kernel t cps =+ let steps g = [kernel t cp g | cp <- cps]+ combine comb step = (\iox -> iox >>= comb) . step+ cycled g = foldl combine return (steps g)+ in cycled++-- | Execute the action once, given the current state in the Markov chain.+execute :: Monad m => Action x m a b -> x -> m (Action x m a b)+execute (viewAction -> Action act fin a) x = liftM (makeAction act fin) (act x a)++-- | Execute the action once every @n@ times, where @n@ is the first argument.+every :: Monad m => Int -> Action x m a b -> Action x m (a,Int) b+every n (viewAction -> Action act fin a) = + let skip_act x (b,i) = if i == (n-1)+ then do b' <- act x b+ return (b',0)+ else return (b,i+1)+ skip_fin = fin . fst+ in makeAction skip_act skip_fin (a,0)
+ MCMC/Distributions.hs view
@@ -0,0 +1,215 @@+{-# LANGUAGE MultiParamTypeClasses, KindSignatures, FlexibleInstances, GADTs, + OverlappingInstances, ViewPatterns #-}++module MCMC.Distributions (+ -- * Methods over distributions+ HasDensity(..)+ , productDensity+ , sampleFrom+ , fromProposal+ -- * Standard distributions+ , uniform+ , mvUniform+ , diag+ , normal+ , mvNormal+ , categorical+ , normalize+ , categoricalNormed+ , beta+ , bern+ , poisson+ ) where++import qualified System.Random.MWC as MWC+import qualified System.Random.MWC.Distributions as MWC.D+import Control.Monad+import qualified Data.Packed.Matrix as M+import qualified Numeric.LinearAlgebra.Algorithms as LA+import qualified Numeric.Container as C+import qualified Language.Hakaru.Distribution as HD+import qualified Language.Hakaru.Types as HT+import Data.Maybe+import Data.Ord+import Data.List as L++import MCMC.Types+import MCMC.SemanticEditors++-- | The class of distributions for which there exists a density method.+class HasDensity d a where+ -- | A method that provides the probability density at a point in the distribution.+ density :: d a -> Density a++-- | Compute the product density of the input distributions.+-- An example use is in constructing a target distribution+-- whose density can be expressed as a product of probability densities.+productDensity :: HasDensity d a => [d a] -> Density a+productDensity ds x = product $ map (\d -> density d x) ds++-- | Get the probability density at a point for any target distribution.+instance HasDensity Target a where+ density (viewTarget -> Target d) = d++-- | Get the probability density at a point for any proposal distribution.+instance HasDensity Proposal a where+ density (viewProposal -> Proposal d _) = d++-- | This function can be used to call the sampling method of any proposal distribution.+sampleFrom :: Proposal a -> Sample a+sampleFrom (viewProposal -> Proposal _ s) = s++{-| + Convenience function for constructing target distributions from predefined+ or custom proposal distributions. One use case is in testing that the samplers+ in the library correctly simulate standard distributions.+ -}+fromProposal :: Proposal a -> Target a+fromProposal = makeTarget . density++-- Uniform -- ++-- Univariate++-- | Univariate uniform distribution over real numbers. The parameters should +-- not be equal. +uniform :: (MWC.Variate a, Real a) => a -> a -> Proposal a+uniform a b+ | b < a = makeUniform b a+ | a < b = makeUniform a b+ | otherwise = error "Wrong parameters for Uniform distribution"++unif1D :: Real a => a -> a -> a -> Double+unif1D a b x+ | x < a = 0+ | x > b = 0+ | otherwise = 1 / realToFrac (b - a)++makeUniform :: (MWC.Variate a, Real a) => a -> a -> Proposal a+makeUniform a b = makeProposal (unif1D a b) (MWC.uniformR (a,b))++-- Multivariate++-- | Multivariate uniform distribution over real numbers.+mvUniform :: (MWC.Variate a, Real a) => [a] -> [a] -> Proposal [a]+mvUniform a b+ | b < a = makeMVUniform b a+ | a < b = makeMVUniform a b+ | otherwise = error "Wrong parameters for multi-variate Uniform distribution"++makeMVUniform :: (MWC.Variate a, Real a) => [a] -> [a] -> Proposal [a]+makeMVUniform a b = + let tupF f (p,q,r) = f p q r+ uniD x = product . map (tupF unif1D) $ zip3 a b x+ uniSF g = mapM (flip MWC.uniformR g) $ zip a b+ in makeProposal uniD uniSF+ +-- Normal --++-- Univariate++-- | Univariate Gaussian distribution over real numbers.+normal :: Double -> Double -> Proposal Double+normal mean cov = makeProposal dens sf+ where hakaruNormal = HD.normal mean (sqrt cov)+ dens x = exp $ HT.logDensity hakaruNormal (HT.Lebesgue x)+ sf g = liftM HT.fromLebesgue $ HT.distSample hakaruNormal g++-- Multivariate++type CovMatrix = M.Matrix Double+type Mu a = M.Matrix a++mu :: M.Element a => [a] -> Mu a+mu mean = M.fromLists [mean]++{-| + Convenience function to create a diagonal matrix from a list representing+ the diagonal. Useful for creating a diagonal covariance matrix for the+ multivariate Gaussian distribution.+ -}+diag :: [Double] -> [[Double]]+diag d = [(nth i) (swapWith e) (replicate (length d) 0) | (i,e) <- zip [1..] d]++-- | Multivariate Gaussian distribution over real numbers.+mvNormal :: [Double] -> [[Double]] -> Proposal [Double]+mvNormal mean cov =+ let covMat = M.fromLists cov+ (muMat, n) = (mu mean, length mean)+ in makeProposal (mvNormalDensity muMat covMat) (mvNormalSF muMat covMat n)++mvNormalDensity :: Mu Double -> CovMatrix -> Density [Double]+mvNormalDensity m cov x = c * exp (-d / 2)+ where (covInv, (lndet, sign)) = LA.invlndet cov+ c1 = (2*pi) ^^ (length x)+ c = 1 / (sqrt $ sign * (exp lndet) * c1)+ xm = C.sub (M.fromLists [x]) m+ prod = xm C.<> covInv C.<> (M.trans xm)+ d = (M.@@>) prod (0,0)++mvNormalSF :: Mu Double -> CovMatrix -> Int -> Sample [Double]+mvNormalSF m cov n g = do+ z <- replicateM n (MWC.D.standard g)+ let zt = M.trans $ M.fromLists [z]+ a = LA.chol cov+ return . head . M.toLists $ C.add m $ C.trans $ a C.<> zt++-- Categorical --++-- | Categorical distribution over instances of the Eq typeclass.+-- The input argument is a list of category-proportion pairs. +-- +-- The input proportions represent relative weights and are not +-- required to be normalized.++-- Look at 'MCMC.Combinators.mixProposals' and 'MCMC.Combinators.mixSteps' for +-- examples of using categorical over elements outside of the 'Eq' typeclass.+categorical :: Eq a => [(a, Double)] -> Proposal a+categorical = categoricalNormed . normalize++-- | Normalize the weights in a list of category-weight pairs.+normalize :: [(a, Double)] -> [(a, Double)]+normalize catProbs = map norm catProbs+ where norm = second (flip (/) s)+ s = foldl (+) 0 $ (snd.unzip) catProbs++-- | Assume that the weights are already normalized. This is useful+-- as an optimized version of @categorical@.+categoricalNormed :: Eq a => [(a,Double)] -> Proposal a+categoricalNormed catProbs =+ -- CHECK: Should fromMaybe default to "error" instead of 0?+ let dens a = snd $ fromMaybe (a,0) $ find ((==)a.fst) catProbs+ in makeProposal dens (categoricalSF catProbs)++categoricalSF :: [(a, Double)] -> Sample a+categoricalSF catProbs g = do+ u <- sampleFrom (uniform 0 1) g+ let (cats, probs) = unzip catProbs+ catsCDF = zip cats $ init $ scanl (+) 0 probs+ return $ fst $ maximumBy (comparing snd) $ filter ((u>=).snd) catsCDF++-- Beta --++-- | Beta distribution over real numbers. Requires non-negative arguments.+beta :: Double -> Double -> Proposal Double+beta a b = makeProposal dens betaSF+ where hakaruBeta = HD.beta a b+ dens x = exp $ HT.logDensity hakaruBeta (HT.Lebesgue x)+ betaSF g = liftM HT.fromLebesgue $ HT.distSample hakaruBeta g++-- Bern --++-- | Bernoulli distribution+bern :: Double -> Proposal Bool+bern p = makeProposal dens bernSF+ where hakaruBern = HD.bern p+ dens x = exp $ HT.logDensity hakaruBern (HT.Discrete x)+ bernSF g = liftM HT.fromDiscrete $ HT.distSample hakaruBern g++-- Poisson --++-- | Univariate Poisson distribution. Requires non-negative argument.+poisson :: Double -> Proposal Int+poisson lambda = makeProposal (exp . HT.logDensity d . HT.Discrete)+ (fmap HT.fromDiscrete . HT.distSample d)+ where d = HD.poisson lambda
+ MCMC/Examples/GMM.hs view
@@ -0,0 +1,475 @@+-- | Sampler for Gaussian Mixture Model+-- +-- Here is the code in the Hakaru language for generating +-- the data used in this example:+-- +-- > p <- unconditioned (beta 2 2)+-- > [m1,m2] <- replicateM 2 $ unconditioned (normal 100 30)+-- > [s1,s2] <- replicateM 2 $ unconditioned (uniform 0 2)+-- > let makePoint = do +-- > b <- unconditioned (bern p)+-- > unconditioned (ifThenElse b (normal m1 s1)+-- > (normal m2 s2))+-- > replicateM nPoints makePoint++module MCMC.Examples.GMM ( GaussianMixtureState (..)+ -- * Target+ -- ** Focus combinators+ -- $focuscombs+ -- , focusLabels, focusGaussParams, focusBernParam, focusObs++ -- ** Record field targets+ -- $fieldtargets+ -- , labelsTarget, gaussParamsTarget, bernParamTarget, obsTarget++ -- ** Target density factors+ -- $targetfactors+ -- , labelsFactor, gaussParamsFactor, bernParamFactor, obsFactor++ -- ** Target density+ -- $tdens+ -- , gmmTarget++ -- * Proposal+ -- ** Proposal update boilerplate+ -- $proposalfocus+ -- , updateLabels, updateGaussParams, updateBernParam++ -- ** Field proposals+ -- $fieldproposals+ -- , labelsProposal, gaussParamsProposal, bernParamProposal++ -- ** Field updaters+ -- $fieldupdaters+ -- , labelsUpdater, gaussParamsUpdater, bernParamUpdater++ -- ** The combined proposal+ -- $gmmprop+ -- , gmmProposal++ -- * Running the sampler+ -- ** Transition kernel+ -- $kernel+ -- , gmmMH++ -- ** Visualization methods+ -- $visual+ -- , histogram, printFields, printLabelN, compareLabels, printHist, batchHist++ -- ** Main+ -- $main+ -- , nPoints, sampleData, gmmStart, gmmTest+ ) where++import MCMC.Combinators+import MCMC.Distributions+import MCMC.Kernels+import MCMC.Actions+import MCMC.Types+import qualified System.Random.MWC as MWC+import qualified Data.Map.Strict as Map+import Control.Monad++data GaussianMixtureState = GMM { labels :: [Bool] -- ^ The list of observation labels+ , gaussParams :: ((Double, Double), (Double, Double)) -- ^ The parameters of the two Gaussians (mean, covariance)+ , bernParam :: Double -- ^ The mixture proportion+ , obs :: [Double] -- ^ The observed data+ }++-- Target ----------++-- Focus combinators++focusLabels :: Target (Double, [Bool]) -> Target GaussianMixtureState+focusLabels t = makeTarget dens+ where dens (GMM l _ p _) = density t (p,l)++focusGaussParams :: Target ((Double, Double), (Double, Double)) -> Target GaussianMixtureState+focusGaussParams t = makeTarget (density t . gaussParams)++focusBernParam :: Target Double -> Target GaussianMixtureState+focusBernParam t = makeTarget (density t . bernParam)++focusObs :: Target ([Bool], ((Double, Double), (Double, Double)), [Double])+ -> Target GaussianMixtureState+focusObs t = makeTarget dens+ where dens (GMM l gps _ o) = density t (l, gps, o)++-- Record field targets++labelsTarget :: Target (Double, [Bool])+labelsTarget = makeTarget $ \(p,ls) -> product $ map (density $ bern p) ls++gaussParamsTarget :: Target ((Double, Double), (Double, Double))+gaussParamsTarget = makeTarget dens+ where dens ((m1, c1), (m2, c2)) = mdens m1 * mdens m2 * cdens c1 * cdens c2+ mdens m = density (normal 100 900) m+ cdens c = density (uniform 0 200) c++bernParamTarget :: Target Double+bernParamTarget = fromProposal (beta 2 2)++obsTarget :: Target ([Bool], ((Double, Double), (Double, Double)), [Double])+obsTarget = makeTarget dens+ where dens (ls, ((m1, c1), (m2, c2)), os) + = let ols = zip os ls+ gauss l = if l then normal m1 (c1*c1) else normal m2 (c2*c2)+ in product $ map (\(o,l) -> density (gauss l) o) ols++-- Target density factors++labelsFactor :: Target GaussianMixtureState+labelsFactor = focusLabels labelsTarget++gaussParamsFactor :: Target GaussianMixtureState+gaussParamsFactor = focusGaussParams gaussParamsTarget++bernParamFactor :: Target GaussianMixtureState+bernParamFactor = focusBernParam bernParamTarget++obsFactor :: Target GaussianMixtureState+obsFactor = focusObs obsTarget++-- Target density++gmmTarget :: Target GaussianMixtureState+gmmTarget = makeTarget $ productDensity + [labelsFactor, gaussParamsFactor, bernParamFactor, obsFactor]++nPoints :: Int+nPoints = 6++sampleData :: [Double]+sampleData = [ 63.13941114139962, 132.02763712240528+ , 62.59642260289356, 132.2616834236893+ , 64.10610391933461, 62.143820541377934 ]++gmmTargetDensityTest :: IO ()+gmmTargetDensityTest = do+ let sampleParams = ((63, 100), (132, 100))+ b = 0.5+ makeState sampleLabels = GMM sampleLabels sampleParams b sampleData+ labels1 = replicate nPoints False+ labels2 = map not labels1+ labels3 = [True, False, True, False, True, True]+ labels4 = [True, True, True, False, False, False]+ putStr $ show labels1 ++ " : " + print $ density gmmTarget $ makeState labels1+ putStr $ show labels2 ++ " : " + print $ density gmmTarget $ makeState labels2+ putStr $ show labels3 ++ " : " + print $ density gmmTarget $ makeState labels3+ putStr $ show labels4 ++ " : " + print $ density gmmTarget $ makeState labels4++-- Proposal ----------++-- Field update combinators++updateLabels :: ([Bool] -> Proposal [Bool]) -> GaussianMixtureState -> Proposal GaussianMixtureState+updateLabels f x = makeProposal dens sf+ where dens y = density (f $ labels x) (labels y)+ sf g = do newLabels <- sampleFrom (f $ labels x) g+ return x { labels = newLabels }++updateGaussParams :: (((Double, Double), (Double, Double)) -> Proposal ((Double, Double), (Double, Double)))+ -> GaussianMixtureState -> Proposal GaussianMixtureState+updateGaussParams f x = makeProposal dens sf+ where dens y = density (f $ gaussParams x) (gaussParams y)+ sf g = do newParams <- sampleFrom (f $ gaussParams x) g+ return x { gaussParams = newParams }++updateBernParam :: (Double -> Proposal Double) -> GaussianMixtureState -> Proposal GaussianMixtureState+updateBernParam f x = makeProposal dens sf+ where dens y = density (f $ bernParam x) (bernParam y)+ sf g = do newParam <- sampleFrom (f $ bernParam x) g+ return x { bernParam = newParam }++-- Field proposals++labelsProposal :: [Bool] -> Proposal [Bool]+labelsProposal ls = chooseProposal nPoints (\n -> updateNth n flipBool ls)+ where flipBool bn = if bn then bern 0 else bern 1++gaussParamsProposal :: ((Double, Double), (Double, Double)) -> Proposal ((Double, Double), (Double, Double))+gaussParamsProposal params = mixProposals $ zip [m1p, c1p, m2p, c2p] (repeat 1)+ where condProp c = normal c 1+ m1p = updateFirst (updateFirst condProp) params+ c1p = updateFirst (updateSecond condProp) params+ m2p = updateSecond (updateFirst condProp) params+ c2p = updateSecond (updateSecond condProp) params++bernParamProposal :: Double -> Proposal Double+bernParamProposal p = uniform (p/2) (1-p/2)++-- Field updaters++labelsUpdater :: GaussianMixtureState -> Proposal GaussianMixtureState+labelsUpdater = updateLabels labelsProposal++gaussParamsUpdater :: GaussianMixtureState -> Proposal GaussianMixtureState+gaussParamsUpdater = updateGaussParams gaussParamsProposal++bernParamUpdater :: GaussianMixtureState -> Proposal GaussianMixtureState+bernParamUpdater = updateBernParam bernParamProposal++-- GMM Proposal++gmmProposal :: GaussianMixtureState -> Proposal GaussianMixtureState+gmmProposal = mixCondProposals $ zip [labelsUpdater, gaussParamsUpdater, bernParamUpdater] [10,1,2]++-- Histogram and other visualizations++histogram :: Ord a => [a] -> Map.Map a Int+histogram ls = foldl addElem Map.empty ls+ where addElem m e = Map.insertWith (+) e 1 m++printFields :: PrintF GaussianMixtureState ([Bool], ((Double, Double), (Double, Double)), Double)+printFields = let f s = (labels s, gaussParams s, bernParam s) in map f ++printLabelN :: Int -> PrintF GaussianMixtureState Bool+printLabelN n = let f s = labels s !! (n-1) in map f++compareLabels :: Int -> Int -> PrintF GaussianMixtureState (Bool,Bool)+compareLabels n m = let f s = (labels s !! (n-1) , labels s !! (m-1)) in map f++printHist :: (Ord s, Show s) => PrintF x s -> Batch x -> IO ()+printHist f (ls,_) = unless (null ls) $ print . histogram $ f ls++batchHist :: (Ord s, Show s) => PrintF x s -> Int -> BatchAction x IO ()+batchHist f n = pack (printHist f) $ inBatches (printHist f) n++-- Kernel ----------++gmmMH :: Step GaussianMixtureState+gmmMH = metropolisHastings gmmTarget gmmProposal++gmmStart :: GaussianMixtureState+gmmStart = GMM { labels = [True, True, True, False, False, False],+ gaussParams = ((63, 100), (132, 100)),+ bernParam = 0.5,+ obs = sampleData }++gmmTest :: IO ()+gmmTest = do+ g <- MWC.createSystemRandom+ let a = batchHist (compareLabels 5 6) 50+ e = every 50 a+ c = every 50 collect+ ls <- walk gmmMH gmmStart (10^6) g c+ putStrLn "Done"+ print $ take 20 (map labels ls)++-- Older, simpler way of writing target density ----------++probLabels :: GaussianMixtureState -> Double+probLabels (GMM l _ p _) = product $ map (\b -> if b then p else 1-p) l++probObs :: GaussianMixtureState -> Double+probObs state = product $ map (\(o,l) -> density (gauss l) o) ols+ where ols = zip (obs state) (labels state)+ ((m1, c1), (m2, c2)) = gaussParams state+ gauss l = if l then normal m1 c1 else normal m2 c2++probGaussParams :: GaussianMixtureState -> Double+probGaussParams state = mdens m1 * mdens m2 * cdens c1 * cdens c2+ where ((m1, c1), (m2, c2)) = gaussParams state+ mdens m = density (normal 100 900) m+ cdens c = density (uniform 0 2) c++probBernParam :: GaussianMixtureState -> Double+probBernParam state = density (beta 2 2) (bernParam state)++gmmTargetOld :: Target GaussianMixtureState+gmmTargetOld = makeTarget dens+ where dens s = probLabels s * probObs s * probGaussParams s * probBernParam s+++-----------------+-- Documentation+-----------------++-- $focuscombs+-- @+-- focusLabels :: Target (Double, [Bool]) -> Target GaussianMixtureState+-- focusLabels t = 'makeTarget' dens+-- where dens (GMM l _ p _) = 'density' t (p,l)+--+-- focusGaussParams :: Target ((Double, Double), (Double, Double)) -> Target GaussianMixtureState+-- focusGaussParams t = 'makeTarget' ('density' t . gaussParams)+--+-- focusBernParam :: Target Double -> Target GaussianMixtureState+-- focusBernParam t = 'makeTarget' ('density' t . bernParam)+--+-- focusObs :: Target ([Bool], ((Double, Double), (Double, Double)), [Double])+-- -> Target GaussianMixtureState+-- focusObs t = 'makeTarget' dens+-- where dens (GMM l gps _ o) = 'density' t (l, gps, o)+-- @+++-- $fieldtargets+-- @+-- labelsTarget :: Target (Double, [Bool])+-- labelsTarget = 'makeTarget' $ \(p,ls) -> product $ map ('density' $ 'bern' p) ls+--+-- gaussParamsTarget :: Target ((Double, Double), (Double, Double))+-- gaussParamsTarget = 'makeTarget' dens+-- where dens ((m1, c1), (m2, c2)) = mdens m1 * mdens m2 * cdens c1 * cdens c2+-- mdens m = 'density' ('normal' 100 900) m+-- cdens c = 'density' ('uniform' 0 200) c+--+-- bernParamTarget :: Target Double+-- bernParamTarget = 'fromProposal' ('beta' 2 2)+--+-- obsTarget :: Target ([Bool], ((Double, Double), (Double, Double)), [Double])+-- obsTarget = 'makeTarget' dens+-- where dens (ls, ((m1, c1), (m2, c2)), os) +-- = let ols = zip os ls+-- gauss l = if l then 'normal' m1 (c1*c1) else 'normal' m2 (c2*c2)+-- in product $ map (\(o,l) -> 'density' (gauss l) o) ols+-- @+++-- $targetfactors+-- @+-- labelsFactor :: Target GaussianMixtureState+-- labelsFactor = focusLabels labelsTarget+--+-- gaussParamsFactor :: Target GaussianMixtureState+-- gaussParamsFactor = focusGaussParams gaussParamsTarget+--+-- bernParamFactor :: Target GaussianMixtureState+-- bernParamFactor = focusBernParam bernParamTarget+--+-- obsFactor :: Target GaussianMixtureState+-- obsFactor = focusObs obsTarget+-- @+++-- $tdens+-- @+-- gmmTarget :: Target GaussianMixtureState+-- gmmTarget = 'makeTarget' $ 'productDensity' +-- [labelsFactor, gaussParamsFactor, bernParamFactor, obsFactor]+-- @+++-- $proposalfocus+-- @+-- updateLabels :: ([Bool] -> Proposal [Bool]) -> GaussianMixtureState -> Proposal GaussianMixtureState+-- updateLabels f x = 'makeProposal' dens sf+-- where dens y = 'density' (f $ labels x) (labels y)+-- sf g = do newLabels <- 'sampleFrom' (f $ labels x) g+-- return x { labels = newLabels }+--+-- updateGaussParams :: (((Double, Double), (Double, Double)) -> Proposal ((Double, Double), (Double, Double)))+-- -> GaussianMixtureState -> Proposal GaussianMixtureState+-- updateGaussParams f x = 'makeProposal' dens sf+-- where dens y = 'density' (f $ gaussParams x) (gaussParams y)+-- sf g = do newParams <- 'sampleFrom' (f $ gaussParams x) g+-- return x { gaussParams = newParams }+--+-- updateBernParam :: (Double -> Proposal Double) -> GaussianMixtureState -> Proposal GaussianMixtureState+-- updateBernParam f x = 'makeProposal' dens sf+-- where dens y = 'density' (f $ bernParam x) (bernParam y)+-- sf g = do newParam <- 'sampleFrom' (f $ bernParam x) g+-- return x { bernParam = newParam }+-- @+++-- $fieldproposals+-- @+-- labelsProposal :: [Bool] -> Proposal [Bool]+-- labelsProposal ls = 'chooseProposal' nPoints (\n -> 'updateNth' n flipBool ls)+-- where flipBool bn = if bn then 'bern' 0 else 'bern' 1+--+-- gaussParamsProposal :: ((Double, Double), (Double, Double)) -> Proposal ((Double, Double), (Double, Double))+-- gaussParamsProposal params = 'mixProposals' $ zip [m1p, c1p, m2p, c2p] (repeat 1)+-- where condProp c = 'normal' c 1+-- m1p = 'updateFirst' ('updateFirst' condProp) params+-- c1p = 'updateFirst' ('updateSecond' condProp) params+-- m2p = 'updateSecond' ('updateFirst' condProp) params+-- c2p = 'updateSecond' ('updateSecond' condProp) params+--+-- bernParamProposal :: Double -> Proposal Double+-- bernParamProposal p = 'uniform' (p/2) (1-p/2)+-- @+++-- $fieldupdaters+-- @+-- labelsUpdater :: GaussianMixtureState -> Proposal GaussianMixtureState+-- labelsUpdater = updateLabels labelsProposal+--+-- gaussParamsUpdater :: GaussianMixtureState -> Proposal GaussianMixtureState+-- gaussParamsUpdater = updateGaussParams gaussParamsProposal+--+-- bernParamUpdater :: GaussianMixtureState -> Proposal GaussianMixtureState+-- bernParamUpdater = updateBernParam bernParamProposal+-- @+++-- $gmmprop+-- @+-- gmmProposal :: GaussianMixtureState -> Proposal GaussianMixtureState+-- gmmProposal = 'mixCondProposals' $ zip [labelsUpdater, gaussParamsUpdater, bernParamUpdater] [10,1,2]+-- @+++-- $kernel+-- @+-- gmmMH :: Step GaussianMixtureState+-- gmmMH = 'metropolisHastings' gmmTarget gmmProposal+-- @+++-- $visual+-- @+-- histogram :: Ord a => [a] -> Map.Map a Int+-- histogram ls = foldl addElem Map.empty ls+-- where addElem m e = Map.insertWith (+) e 1 m+--+-- printFields :: PrintF GaussianMixtureState ([Bool], ((Double, Double), (Double, Double)), Double)+-- printFields = let f s = (labels s, gaussParams s, bernParam s) in map f +--+-- printLabelN :: Int -> PrintF GaussianMixtureState Bool+-- printLabelN n = let f s = labels s !! (n-1) in map f+--+-- compareLabels :: Int -> Int -> PrintF GaussianMixtureState (Bool,Bool)+-- compareLabels n m = let f s = (labels s !! (n-1) , labels s !! (m-1)) in map f+--+-- printHist :: (Ord s, Show s) => PrintF x s -> Batch x -> IO ()+-- printHist f (ls,_) = unless (null ls) $ print . histogram $ f ls+--+-- batchHist :: (Ord s, Show s) => PrintF x s -> Int -> BatchAction x IO ()+-- batchHist f n = 'pack' (printHist f) $ 'inBatches' (printHist f) n+-- @+++-- $main+-- @+-- nPoints :: Int+-- nPoints = 6+--+-- sampleData :: [Double]+-- sampleData = [ 63.13941114139962, 132.02763712240528+-- , 62.59642260289356, 132.2616834236893+-- , 64.10610391933461, 62.143820541377934 ]+--+-- gmmStart :: GaussianMixtureState+-- gmmStart = GMM { labels = [True, True, True, False, False, False],+-- gaussParams = ((63, 100), (132, 100)),+-- bernParam = 0.5,+-- obs = sampleData }+--+-- gmmTest :: IO ()+-- gmmTest = do+-- g <- MWC.createSystemRandom+-- let a = batchHist (compareLabels 5 6) 50+-- e = 'every' 50 a+-- c = 'every' 50 'collect'+-- ls <- 'walk' gmmMH gmmStart (10^6) g c+-- putStrLn \"Done\"+-- print $ take 20 (map labels ls)+-- @
+ MCMC/Examples/HandwrittenGMM.hs view
@@ -0,0 +1,311 @@+-- | Optimized sampler for Gaussian Mixture Model+-- +-- Here is the code in the Hakaru language for generating +-- the data used in this example:+-- +-- > p <- unconditioned (beta 2 2)+-- > [m1,m2] <- replicateM 2 $ unconditioned (normal 100 30)+-- > [s1,s2] <- replicateM 2 $ unconditioned (uniform 0 2)+-- > let makePoint = do +-- > b <- unconditioned (bern p)+-- > unconditioned (ifThenElse b (normal m1 s1)+-- > (normal m2 s2))+-- > replicateM nPoints makePoint++module MCMC.Examples.HandwrittenGMM (+ GaussianMixtureState(..)+ -- * Focused targets+ -- $tar++ -- * Focused proposals+ -- $prop+ + -- * Focused steps+ -- *** Each step computes only those parts of the density ratio that its proposal affects - the other parts would cancel out+ -- $steps+ + -- * Optimized sampler+ -- *** A mixture of focused, i.e. optimized steps + -- $sampler+ + -- * Main+ -- $main+ ) where++import MCMC.Types+import MCMC.Kernels+import MCMC.Distributions+import MCMC.Actions+import MCMC.Combinators+import qualified System.Random.MWC as MWC++data GaussianMixtureState = GMM { labels :: [Bool]+ , gaussParams :: ((Double, Double), (Double, Double)) + , bernParam :: Double }++nPoints :: Int+nPoints = 6++stepLabels :: [Double] -> Step GaussianMixtureState+stepLabels obs = chooseStep nPoints + (\i -> makeTarget $ dens i) labelsProposal metropolisHastings+ where dens i state = density (targetLabel i) state *+ density (targetObs i obs) state++-- This could be optimized further if we know the label corresponding +-- to the gaussian to which the updated param belongs.+stepGaussParams :: [Double] -> Step GaussianMixtureState+stepGaussParams obs = metropolisHastings (makeTarget dens) gaussParamsProposal+ where dens state = density targetGaussParams state * + product [density (targetObs i obs) state | i <- [1..nPoints]]++stepBernParam :: Step GaussianMixtureState+stepBernParam = metropolisHastings (makeTarget dens) bernParamProposal+ where dens state = density targetBernParam state *+ product [density (targetLabel i) state | i <- [1..nPoints]]++gmmSampler :: [Double] -> Step GaussianMixtureState+gmmSampler obs = mixSteps $ + zip [(stepLabels obs), (stepGaussParams obs), stepBernParam] [1,1,1] ++-- | Main++sampleData :: [Double]+sampleData = [ 63.13941114139962, 132.02763712240528+ , 62.59642260289356, 132.2616834236893+ , 64.10610391933461, 62.143820541377934 ]++startState :: GaussianMixtureState+startState = GMM { -- labels = [True, True, True, False, False, False],+ labels = [False, False, False, True, True, True],+ gaussParams = ((63, 100), (132, 100)),+ bernParam = 0.5 }++test :: IO ()+test = do+ g <- MWC.createSystemRandom+ let c = every 50 collect+ p = every 1 (display labels)+ -- ls <- walk (gmmSampler sampleData) startState (10^6) g c+ -- print $ take 20 (map labels ls)+ walk (gmmSampler sampleData) startState (10^2) g p+ +-- | Targets -- ++-- Labels++targetLabel :: Int -> Target GaussianMixtureState+targetLabel i = makeTarget (densityLabel i)++densityLabel :: Int -> GaussianMixtureState -> Double+densityLabel i (GMM l _ p) = if (l !! (i-1)) then p else 1-p ++-- Gauss params++targetGaussParams :: Target GaussianMixtureState+targetGaussParams = makeTarget densityGaussParams++densityGaussParams :: GaussianMixtureState -> Double+densityGaussParams state = mdens m1 * mdens m2 * cdens c1 * cdens c2+ where ((m1, c1), (m2, c2)) = gaussParams state+ mdens m = density (normal 100 900) m+ cdens c = density (uniform 0 2) c++-- Bern param++targetBernParam :: Target GaussianMixtureState+targetBernParam = makeTarget densityBernParam++densityBernParam :: GaussianMixtureState -> Double+densityBernParam state = density (beta 2 2) (bernParam state)++-- Obs / data points++targetObs :: Int -> [Double] -> Target GaussianMixtureState+targetObs i obs = makeTarget (densityObs i obs)++densityObs :: Int -> [Double] -> GaussianMixtureState -> Double+densityObs i obs state = if labels state !! (i-1)+ then density (normal m1 c1) oi+ else density (normal m2 c2) oi+ where oi = obs !! (i-1)+ ((m1, c1), (m2, c2)) = gaussParams state++-- | Proposals --++-- Labels++labelsProposal :: Int -> GaussianMixtureState -> Proposal GaussianMixtureState+labelsProposal i x = makeProposal dens sf+ where dens y = density (updateLabel i $ labels x) (labels y)+ sf g = do newLabels <- sampleFrom (updateLabel i $ labels x) g+ return x { labels = newLabels }++updateLabel :: Int -> [Bool] -> Proposal [Bool]+updateLabel i ls = updateNth i flipBool ls+ where flipBool bn = if bn then bern 0 else bern 1++-- Gauss params++gaussParamsProposal :: GaussianMixtureState -> Proposal GaussianMixtureState+gaussParamsProposal x = makeProposal dens sf+ where dens y = density (updateGaussParams $ gaussParams x) (gaussParams y)+ sf g = do newParams <- sampleFrom (updateGaussParams $ gaussParams x) g+ return x { gaussParams = newParams }++updateGaussParams :: ((Double, Double), (Double, Double)) -> Proposal ((Double, Double), (Double, Double))+updateGaussParams params = mixProposals $ zip [m1p, c1p, m2p, c2p] (repeat 1)+ where condProp c = normal c 1+ m1p = updateFirst (updateFirst condProp) params+ c1p = updateFirst (updateSecond condProp) params+ m2p = updateSecond (updateFirst condProp) params+ c2p = updateSecond (updateSecond condProp) params++-- Bern param++bernParamProposal :: GaussianMixtureState -> Proposal GaussianMixtureState+bernParamProposal x = makeProposal dens sf+ where dens y = density (updateBernParam $ bernParam x) (bernParam y)+ sf g = do newParam <- sampleFrom (updateBernParam $ bernParam x) g+ return x { bernParam = newParam }++updateBernParam :: Double -> Proposal Double+updateBernParam p = uniform (p/2) (1-p/2)+++-----------------+-- Documentation+-----------------++-- $tar+-- @+-- targetLabel :: Int -> Target GaussianMixtureState+-- targetLabel i = 'makeTarget' (densityLabel i)+--+-- densityLabel :: Int -> GaussianMixtureState -> Double+-- densityLabel i (GMM l _ p) = if (l !! (i-1)) then p else 1-p +--+--+-- targetGaussParams :: Target GaussianMixtureState+-- targetGaussParams = 'makeTarget' densityGaussParams+--+-- densityGaussParams :: GaussianMixtureState -> Double+-- densityGaussParams state = mdens m1 * mdens m2 * cdens c1 * cdens c2+-- where ((m1, c1), (m2, c2)) = gaussParams state+-- mdens m = 'density' ('normal' 100 900) m+-- cdens c = 'density' ('uniform' 0 2) c+--+-- +-- targetBernParam :: Target GaussianMixtureState+-- targetBernParam = 'makeTarget' densityBernParam+--+-- densityBernParam :: GaussianMixtureState -> Double+-- densityBernParam state = 'density' ('beta' 2 2) (bernParam state)+--+--+-- targetObs :: Int -> [Double] -> Target GaussianMixtureState+-- targetObs i obs = 'makeTarget' (densityObs i obs)+--+-- densityObs :: Int -> [Double] -> GaussianMixtureState -> Double+-- densityObs i obs state = if labels state !! (i-1)+-- then 'density' ('normal' m1 c1) oi+-- else 'density' ('normal' m2 c2) oi+-- where oi = obs !! (i-1)+-- ((m1, c1), (m2, c2)) = gaussParams state+-- @+++-- $prop+-- @+-- labelsProposal :: Int -> GaussianMixtureState -> Proposal GaussianMixtureState+-- labelsProposal i x = 'makeProposal' dens sf+-- where dens y = 'density' (updateLabel i $ labels x) (labels y)+-- sf g = do newLabels <- 'sampleFrom' (updateLabel i $ labels x) g+-- return x { labels = newLabels }+--+-- updateLabel :: Int -> [Bool] -> Proposal [Bool]+-- updateLabel i ls = 'updateNth' i flipBool ls+-- where flipBool bn = if bn then 'bern' 0 else 'bern' 1+--+--+-- gaussParamsProposal :: GaussianMixtureState -> Proposal GaussianMixtureState+-- gaussParamsProposal x = 'makeProposal' dens sf+-- where dens y = 'density' (updateGaussParams $ gaussParams x) (gaussParams y)+-- sf g = do newParams <- 'sampleFrom' (updateGaussParams $ gaussParams x) g+-- return x { gaussParams = newParams }+--+-- updateGaussParams :: ((Double, Double), (Double, Double)) -> Proposal ((Double, Double), (Double, Double))+-- updateGaussParams params = 'mixProposals' $ zip [m1p, c1p, m2p, c2p] (repeat 1)+-- where condProp c = 'normal' c 1+-- m1p = 'updateFirst' ('updateFirst' condProp) params+-- c1p = 'updateFirst' ('updateSecond' condProp) params+-- m2p = 'updateSecond' ('updateFirst' condProp) params+-- c2p = 'updateSecond' ('updateSecond' condProp) params+--+--+-- bernParamProposal :: GaussianMixtureState -> Proposal GaussianMixtureState+-- bernParamProposal x = 'makeProposal' dens sf+-- where dens y = 'density' (updateBernParam $ bernParam x) (bernParam y)+-- sf g = do newParam <- 'sampleFrom' (updateBernParam $ bernParam x) g+-- return x { bernParam = newParam }+--+-- updateBernParam :: Double -> Proposal Double+-- updateBernParam p = 'uniform' (p/2) (1-p/2)+-- @+++-- $steps+-- @+-- stepLabels :: [Double] -> Step GaussianMixtureState+-- stepLabels obs = 'chooseStep' nPoints +-- (\i -> 'makeTarget' $ dens i) labelsProposal 'metropolisHastings'+-- where dens i state = 'density' (targetLabel i) state *+-- 'density' (targetObs i obs) state+--+-- -- This could be optimized further if we know the label corresponding +-- -- to the gaussian to which the updated param belongs.+-- stepGaussParams :: [Double] -> Step GaussianMixtureState+-- stepGaussParams obs = 'metropolisHastings' ('makeTarget' dens) gaussParamsProposal+-- where dens state = 'density' targetGaussParams state * +-- product ['density' (targetObs i obs) state | i <- [1..nPoints]]+--+-- stepBernParam :: Step GaussianMixtureState+-- stepBernParam = 'metropolisHastings' ('makeTarget' dens) bernParamProposal+-- where dens state = 'density' targetBernParam state *+-- product ['density' (targetLabel i) state | i <- [1..nPoints]]+-- @+++-- $sampler+-- @+-- gmmSampler :: [Double] -> Step GaussianMixtureState+-- gmmSampler obs = 'mixSteps' $ +-- zip [(stepLabels obs), (stepGaussParams obs), stepBernParam] [1,1,1]+-- @+++-- $main+-- @+-- nPoints :: Int+-- nPoints = 6+--+-- sampleData :: [Double]+-- sampleData = [ 63.13941114139962, 132.02763712240528+-- , 62.59642260289356, 132.2616834236893+-- , 64.10610391933461, 62.143820541377934 ]+--+-- startState :: GaussianMixtureState+-- startState = GMM { -- labels = [True, True, True, False, False, False],+-- labels = [False, False, False, True, True, True],+-- gaussParams = ((63, 100), (132, 100)),+-- bernParam = 0.5 }+--+-- test :: IO ()+-- test = do+-- g <- MWC.createSystemRandom+-- let c = 'every' 50 'collect'+-- p = 'every' 1 ('display' labels)+-- -- ls <- 'walk' (gmmSampler sampleData) startState (10^6) g c+-- -- print $ take 20 (map labels ls)+-- 'walk' (gmmSampler sampleData) startState (10^2) g p+-- @
+ MCMC/Kernels.hs view
@@ -0,0 +1,111 @@+{-# LANGUAGE GADTs, MultiParamTypeClasses, KindSignatures, FlexibleInstances,+ ViewPatterns #-}++module MCMC.Kernels ( + -- * Making the random walk+ walk+ -- * Transition kernels+ -- ** Metropolis-Hastings+ , metropolisHastings+ , vizMH+ , printMH+ -- ** Simulated Annealing+ , Temp+ , CoolingSchedule+ , StateSA+ , simulatedAnnealing+ , vizSA+ , printSA+ -- ** Gibbs+ , gibbs+ ) where++import MCMC.Actions+import MCMC.Combinators+import MCMC.Distributions+import MCMC.Types+import Text.Printf++-- | Execute a random walk and create a Markov chain.+walk :: Step x -- ^ The stepping style (based on the transition kernel) + -> x -- ^ The starting state+ -> Int -- ^ The number of steps to take+ -> Rand -- ^ A PRNG+ -> Action x IO a b -- ^ An action to take at each step in the walk+ -> IO b -- ^ The action-dependent output at the end of the walk+walk _ _ 0 _ (viewAction -> Action _ f a) = f a+walk step x n r action = do + x' <- step r x+ execute action x' >>= walk step x' (n-1) r++-- Metropolis Hastings --++metropolisHastings :: Kernel a a+metropolisHastings t c_p = + let mhStep g xi = do+ u <- sampleFrom (uniform 0 1) g+ xstar <- sampleFrom (c_p xi) g+ let accept = min 1 (numer / denom)+ numer = density t xstar * density (c_p xstar) xi+ denom = density t xi * density (c_p xi) xstar+ return $ if u < accept then xstar else xi+ in mhStep++vizMH :: PrintF Double Double+vizMH = id++-- Visualizes only the first dimension+vizMHFirstDim :: PrintF [Double] Double+vizMHFirstDim = map head++printMH :: PrintF [Double] [String]+printMH lls = let p d = printf "%0.3f" d :: String+ in map (map p) lls++-- Simulated Annealing --++type Temp = Double+-- | This is the tempering function used in the simulated annealing process.+type CoolingSchedule = Temp -> Temp+type StateSA a = (a, Temp, CoolingSchedule)++simulatedAnnealing :: Kernel (StateSA a) a+simulatedAnnealing t c_p = + let saStep g (xi,temp,cool) = do+ u <- sampleFrom (uniform 0 1) g+ xstar <- sampleFrom (c_p xi) g+ let accept = min 1 (numer / denom)+ numer = (*) (density (c_p xstar) xi) $ (**) (density t xstar) (1 / temp)+ denom = (*) (density (c_p xi) xstar) $ (**) (density t xi) (1 / temp)+ new_temp = cool temp+ return $ if u < accept then (xstar,new_temp,cool) else (xi,new_temp,cool)+ in saStep++tripleFirst :: (a, b, c) -> a+tripleFirst (a,_,_) = a++myFilter :: [[Double]] -> [[Double]]+myFilter = filter (\x -> x < (repeat 15) && x > (repeat $ -5))++vizSA :: PrintF (StateSA Double) Double+vizSA = map tripleFirst++-- Visualizes only the first dimension+vizSAFirstDim :: PrintF (StateSA [Double]) Double+vizSAFirstDim = vizMHFirstDim . myFilter . map tripleFirst++-- Print the samples without the temp/cooling schedule+printSA :: PrintF (StateSA [Double]) [String]+printSA = printMH . myFilter . map tripleFirst++-- Gibbs --++-- | The full conditional proposals must be specified to this transition kernel.+gibbs :: Target a -> [a -> Proposal a] -> Step a+gibbs = cycleStep alwaysAccept+ where alwaysAccept _ q g x = sampleFrom (q x) g++-- MCMC-EM --++mhEM :: Kernel theta theta+mhEM t c_p = undefined
+ MCMC/SemanticEditors.hs view
@@ -0,0 +1,74 @@+-- | Semantic editor combinators+-- +-- <http://conal.net/blog/posts/semantic-editor-combinator>++module MCMC.SemanticEditors ( first+ , second+ , car+ , cdr+ , nth+ , nthM+ , block+ , blockM+ , swapWith+ , chopAt+ ) where++first :: (a -> a') -> ((a,b) -> (a',b))+second :: (b -> b') -> ((a,b) -> (a,b'))+ +first f = \ (a,b) -> (f a, b)+second g = \ (a,b) -> (a, g b)++car :: (a -> a) -> ([a] -> [a])+cdr :: ([a] -> [a]) -> ([a] -> [a])++car f = \(x:xs) -> f x : xs+cdr f = \(x:xs) -> x : f xs++-- A sample is 1-indexed, i.e., dimensions go from 1 to n++nth :: Int -> (a -> a) -> ([a] -> [a])+nth 1 = car+nth n = cdr . nth (n-1)++carM :: Monad m => (a -> m a) -> ([a] -> m [a])+carM f (x:xs) = do x' <- f x + return $ x' : xs++cdrM :: Monad m => ([a] -> m [a]) -> ([a] -> m [a])+cdrM f (x:xs) = do xs' <- f xs+ return $ x : xs'++nthM :: Monad m => Int -> (a -> m a) -> ([a] -> m [a])+nthM 1 = carM+nthM n = cdrM . nthM (n-1)++block :: Int -> Int -> ([a] -> [a]) -> ([a] -> [a])+block begin end f ls+ | begin == end = nth begin (single f) ls+ | begin == 1 = f (take end ls) ++ drop end ls+ | otherwise = front ++ f mids ++ back+ where (front, mids, back) = chopAt begin end ls++blockM :: Monad m => Int -> Int -> ([a] -> m [a]) -> ([a] -> m [a])+blockM begin end f ls+ | begin == end = nthM begin (singleM f) ls+ | begin == 1 = f (take end ls) >>= return . flip (++) (drop end ls)+ | otherwise = do let (front, mids, rest) = chopAt begin end ls+ fmids <- f mids+ return $ front ++ fmids ++ rest++swapWith :: a -> (b -> a)+swapWith x _ = x++single :: ([a] -> [a]) -> a -> a+single f x = head $ f [x]++singleM :: Monad m => ([a] -> m [a]) -> a -> m a+singleM f x = f [x] >>= return.head++chopAt :: Int -> Int -> [a] -> ([a], [a], [a])+chopAt begin end ls = (front, mids, back)+ where (front, rest) = splitAt (begin-1) ls+ (mids, back) = splitAt (end + 1 - begin) rest
+ MCMC/Tests.hs view
@@ -0,0 +1,125 @@+module MCMC.Tests where++import MCMC.Types+import MCMC.Distributions+import MCMC.Kernels+import MCMC.Actions+import MCMC.Combinators+import qualified System.Random.MWC as MWC++-- Bimodal distribution from section 3.1 of+-- "An Introduction to MCMC for Machine Learning" by C. Andrieu et al.+exampleTarget :: Target Double+exampleTarget = + makeTarget $ \x -> 0.3 * exp (-0.2*x*x) + 0.7 * exp (-0.2 * (x-10)**2)++gaussianProposal :: Double -> Proposal Double+gaussianProposal x = normal x 10000++exampleMH :: Step Double+exampleMH = metropolisHastings exampleTarget gaussianProposal++mhTest :: IO ()+mhTest = do+ g <- MWC.createSystemRandom+ let a = batchViz vizMH 50+ e = every 100 a+ walk exampleMH 0 (10^6) g e++exampleSA :: Step (StateSA Double)+exampleSA = simulatedAnnealing exampleTarget gaussianProposal++saTest :: IO ()+saTest = do+ g <- MWC.createSystemRandom+ let coolSch = (*) (1 - 1e-3) :: Temp -> Temp+ x0 = (0, 1, coolSch)+ a = batchViz vizSA 50+ e = every 100 a+ walk exampleSA x0 (10^6) g e++gMix :: Target [Double]+gMix = let g1 = mvNormal [0,0] (diag [1,1])+ g2 = mvNormal [5,5] (diag [2,2])+ in mixTargets $ zip [g1, g2] [0.3, 0.7] ++prop1 :: [Double] -> Proposal [Double]+prop1 x = updateNth 1 (\y -> normal y 1) x++prop2 :: [Double] -> Proposal [Double]+prop2 x = updateNth 2 (\y -> normal y 1) x++mhMix :: Step [Double]+mhMix = let mh1 = metropolisHastings gMix prop1+ mh2 = metropolisHastings gMix prop2+ in mixSteps $ zip [mh1, mh2] [0.7, 0.3] ++mixTest :: IO ()+mixTest = do+ g <- MWC.createSystemRandom+ let a = batchPrint printMH 50+ walk mhMix [0,0] (10^6) g a++mhCycle :: Step [Double]+mhCycle = cycleStep metropolisHastings gMix [prop1, prop2]++cycleTest :: IO ()+cycleTest = do+ g <- MWC.createSystemRandom+ let a = batchPrint printMH 50+ walk mhCycle [0,0] (10^6) g a++blockMH :: Step [Double]+blockMH = let target = fromProposal $ mvNormal [0,1,4,7] (diag [2,2,2,2])+ mh4D = metropolisHastings target+ mh1 = mh4D $ updateBlock 1 2 (\y -> mvNormal y (diag [1,1]))+ mh2 = mh4D $ updateBlock 3 3 (\y -> mvNormal y [[1]])+ mh3 = mh4D $ updateNth 4 (\y -> normal y 1)+ in mixSteps $ zip [mh1, mh2, mh3] [0.5, 0.4, 0.7] ++blockTest :: IO ()+blockTest = do+ g <- MWC.createSystemRandom+ let a = batchPrint printMH 50+ walk blockMH [0,0,0,0] (10^6) g a++-- main :: IO ()+-- main = blockTest++gTest :: IO ()+gTest = do + print $ density (prop1 [1,2]) [2,4]+ +gaussMix :: Proposal Double+gaussMix = mixProposals $ zip [normal 0 1, normal 5 20, normal 10 0.1] [1..]++betaMix :: Proposal Double+betaMix = mixProposals $ zip [beta 1 1, beta 2 2, beta 3 3] [1..]++mixOfMixes :: Proposal Double+mixOfMixes = mixProposals [(gaussMix, 3), (betaMix, 1)]++lol :: [[Double]] -> Proposal [[Double]]+lol = updateNth 4 (updateNth 2 (\y -> normal y 2))++betaUpdate :: [Double] -> Proposal [Double]+betaUpdate = updateNth 1 (\x -> beta x 3)++gaussUpdate :: [Double] -> Proposal [Double]+gaussUpdate = updateBlock 1 2 (\y -> mvNormal y (diag [1,1]))++mhMix2 :: Step [Double]+mhMix2 = let mh1 = metropolisHastings gMix betaUpdate+ mh2 = metropolisHastings gMix gaussUpdate+ in mixSteps $ zip [mh1, mh2] [0.7, 0.3]++bimodal :: Target [Double]+bimodal = makeTarget dens+ where dens [x,y] = 0.3 * exp (-0.2*x*x) + 0.7 * exp (-0.2 * (y-10)**2)++mhSampler :: Step [Double]+mhSampler = metropolisHastings bimodal betaUpdate++toggle :: (Double, (Bool,String)) -> Proposal (Double, (Bool,String))+toggle = updateSecond (updateFirst (\ b -> bern $ if b then 0 else 1))+
+ MCMC/Types.hs view
@@ -0,0 +1,110 @@+{-# LANGUAGE MultiParamTypeClasses #-}++module MCMC.Types ( + -- * Targets and Proposals+ Rand+ , Density+ , Sample+ , Target+ , viewTarget+ , TargetView(..)+ , makeTarget+ , Proposal+ , viewProposal+ , ProposalView(..)+ , makeProposal+ -- * Transition kernels+ , Step+ , Kernel+ -- * Actions+ , Act+ , Action+ , viewAction+ , ActionView(..)+ , makeAction+ ) where++import qualified System.Random.MWC as MWC++-- | An even shorter name for PRNGs in the 'IO' monad.+type Rand = MWC.GenIO++-- | The probability density function used in both target and proposal distributions.+-- Given an input point, this method returns a probability density. +type Density a = a -> Double++-- | The type for target distributions that can be used in any MCMC sampler.+newtype Target a = T {viewTarget :: TargetView a}+newtype TargetView a = Target (Density a)++-- | Method for constructing custom target distributions.+-- +-- Target distributions need only a density method.+makeTarget :: Density a -> Target a+makeTarget = T . Target++-- | A procedure that, given a source of randomness, returns an action that +-- produces a sample. The type itself is read as a verb, i.e, "to sample".+type Sample a = Rand -> IO a++-- | The type for proposal distributions that can be used in any MCMC +-- sampler.+newtype Proposal a = P {viewProposal :: ProposalView a}+data ProposalView a = Proposal (Density a) (Sample a)++-- | Method for constructing custom proposal distributions.+-- +-- Proposal distributions need both a density and a sampling method.+makeProposal :: Density a -> Sample a -> Proposal a+makeProposal d = P . Proposal d++-- Kernels --++-- | The type for one step in the random walk. +-- A value of type @'Step' a@ is a function that takes a source of randomness and a +-- current state and returns an action producing a subsequent state.+type Step x = Rand -> x -> IO x++-- | The type for MCMC transition kernels. +-- +-- The input arguments are the target +-- distribution (to be modeled) and a /conditional/ proposal distribution.+-- +-- The result is a 'Step' that will make one move in the random walk based+-- on the current state.+-- In general, an MCMC kernel consists of using:+-- +-- * the conditioned proposal to make a hypothesis move, and then+-- * the semantics of the specific 'Kernel' at hand to either accept +-- this hypothesis (and move to the new+-- state), or reject the hypothesis (and stay at the current state).+-- +-- Type parameter definitions:+-- +-- [@x@] The kernel-state. This is the type for each state in the /Markov chain/.+-- [@a@] The distribution-state. This is the domain of the target distribution as+-- well as the type of values sampled from the proposal distribution.+-- +-- In general, we need different types to represent the kernel-state and +-- distribution-state because the+-- kernel-state may hold extra information that gets updated with each step.+-- Look at 'MCMC.Kernels.simulatedAnnealing' for+-- an example where @x@ differs from @a@.+type Kernel x a = Target a -> (a -> Proposal a) -> Step x++type Act x m a = x -> a -> m a++-- | Type parameter definitions:+-- +-- [@x@] The kernel-state (see 'MCMC.Types.Kernel')+-- [@a@] The action-state, specific to the action being performed+-- [@m@] The monad in which the action is performed+-- [@b@] The final returned state type+newtype Action x m a b = A {viewAction :: ActionView x m a b}+data ActionView x m a b = Action (Act x m a) (a -> m b) a++makeAction :: Act x m a -- ^ The action to perform at each step of the random walk+ -> (a -> m b) -- ^ The /finish/ function, called at the end of the sampling process+ -> a -- ^ The current value of the action-state+ -> Action x m a b+makeAction act fin = A . (Action act fin)
− Tests.hs
@@ -1,88 +0,0 @@-module Main where--import Distributions-import Kernels-import Actions-import qualified System.Random.MWC as MWC---- Bimodal distribution from section 3.1 of--- "An Introduction to MCMC for Machine Learning" by C. Andrieu et al.-exampleTarget :: Target [Double]-exampleTarget = - T $ \[x] -> 0.3 * exp (-0.2*x*x) + 0.7 * exp (-0.2 * (x-10)**2)--gaussianProposal :: [Double] -> Proposal [Double]-gaussianProposal x = normal x [[10000]]--exampleMH :: Step [Double]-exampleMH = metropolisHastings exampleTarget gaussianProposal--mhTest :: IO ()-mhTest = do- g <- MWC.createSystemRandom- let a = batchViz vizMH 50- e = every 100 a- walk exampleMH [0] (10^6) g e--exampleSA :: Step (StateSA [Double])-exampleSA = simulatedAnnealing exampleTarget gaussianProposal--saTest :: IO ()-saTest = do- g <- MWC.createSystemRandom- let coolSch = (*) (1 - 1e-3) :: Temp -> Temp- x0 = ([0], 1, coolSch)- a = batchViz vizSA 50- e = every 100 a- walk exampleSA x0 (10^6) g e--gMix :: Target [Double]-gMix = let g1 = normal [0,0] (diag [1,1])- g2 = normal [5,5] (diag [2,2])- in targetMix [g1, g2] [0.3, 0.7] --prop1 :: [Double] -> Proposal [Double]-prop1 x = updateNth 1 (\y -> normal y [[1]]) x--prop2 :: [Double] -> Proposal [Double]-prop2 x = updateNth 2 (\y -> normal y [[1]]) x--mhMix = let mh1 = metropolisHastings gMix prop1- mh2 = metropolisHastings gMix prop2- in mixSteps [mh1, mh2] [0.7, 0.3] --mixTest :: IO ()-mixTest = do- g <- MWC.createSystemRandom- let a = batchPrint printMH 50- walk mhMix [0,0] (10^6) g a--mhCycle :: Step [Double]-mhCycle = cycleKernel metropolisHastings gMix [prop1, prop2]--cycleTest :: IO ()-cycleTest = do- g <- MWC.createSystemRandom- let a = batchPrint printMH 50- walk mhCycle [0,0] (10^6) g a--blockMH :: Step [Double]-blockMH = let target = fromProposal $ normal [0,1,4,7] (diag [2,2,2,2])- mh4D = metropolisHastings target- mh1 = mh4D $ updateBlock 1 2 (\y -> normal y (diag [1,1]))- mh2 = mh4D $ updateBlock 3 3 (\y -> normal y [[1]])- mh3 = mh4D $ updateNth 4 (\y -> normal y [[1]])- in mixSteps [mh1, mh2, mh3] [0.5, 0.4, 0.7] --blockTest :: IO ()-blockTest = do- g <- MWC.createSystemRandom- let a = batchPrint printMH 50- walk blockMH [0,0,0,0] (10^6) g a--main :: IO ()-main = blockTest--gTest :: IO ()-gTest = do - print $ density (prop1 [1,2]) [2,4]
mcmc-samplers.cabal view
@@ -2,23 +2,31 @@ -- documentation, see http://haskell.org/cabal/users-guide/ name: mcmc-samplers-version: 0.1.0.0-synopsis: A library of combinators to build MCMC kernels, proposals, and targets--- description: -license: PublicDomain+version: 0.1.1.0+synopsis: Combinators for MCMC sampling+description: A library of combinators to build transition kernels, proposal distributions, target distributions, and stream operations for MCMC sampling.+license: BSD3 license-file: LICENSE author: Praveen Narayanan maintainer: revenap@gmail.com -- copyright: category: Math build-type: Simple-extra-source-files: Gibbs.hs, Tests.hs, README.md+extra-source-files: Gibbs.hs, README.md cabal-version: >=1.10 library- exposed-modules: Actions, Distributions, Kernels- -- other-modules: + exposed-modules: MCMC.Actions, MCMC.Combinators, MCMC.Distributions,+ MCMC.Kernels, MCMC.Types, MCMC.SemanticEditors, + MCMC.Examples.GMM, MCMC.Examples.HandwrittenGMM+ other-modules: MCMC.Tests other-extensions: MultiParamTypeClasses, KindSignatures, FlexibleInstances, GADTs, OverlappingInstances- build-depends: base >=4.6 && <4.7, mwc-random >=0.13 && <0.14, primitive >=0.5 && <0.6, hmatrix >=0.15 && <0.16, vector >=0.10 && <0.11, statistics >=0.11 && <0.12, containers >=0.5 && <0.6+ build-depends: base >=4.6 && <5,+ mwc-random >=0.13 && <0.14,+ primitive >=0.5 && <0.6,+ hmatrix >=0.15, + statistics >=0.11,+ containers >=0.5 && <0.6,+ hakaru >=0.1.3 -- hs-source-dirs: default-language: Haskell2010