diff --git a/Actions.hs b/Actions.hs deleted file mode 100644 --- a/Actions.hs +++ /dev/null @@ -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 - diff --git a/Distributions.hs b/Distributions.hs deleted file mode 100644 --- a/Distributions.hs +++ /dev/null @@ -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 diff --git a/Kernels.hs b/Kernels.hs deleted file mode 100644 --- a/Kernels.hs +++ /dev/null @@ -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 - diff --git a/MCMC/Actions.hs b/MCMC/Actions.hs new file mode 100644 --- /dev/null +++ b/MCMC/Actions.hs @@ -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 + diff --git a/MCMC/Combinators.hs b/MCMC/Combinators.hs new file mode 100644 --- /dev/null +++ b/MCMC/Combinators.hs @@ -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) diff --git a/MCMC/Distributions.hs b/MCMC/Distributions.hs new file mode 100644 --- /dev/null +++ b/MCMC/Distributions.hs @@ -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 diff --git a/MCMC/Examples/GMM.hs b/MCMC/Examples/GMM.hs new file mode 100644 --- /dev/null +++ b/MCMC/Examples/GMM.hs @@ -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) +-- @ diff --git a/MCMC/Examples/HandwrittenGMM.hs b/MCMC/Examples/HandwrittenGMM.hs new file mode 100644 --- /dev/null +++ b/MCMC/Examples/HandwrittenGMM.hs @@ -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 +-- @ diff --git a/MCMC/Kernels.hs b/MCMC/Kernels.hs new file mode 100644 --- /dev/null +++ b/MCMC/Kernels.hs @@ -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 diff --git a/MCMC/SemanticEditors.hs b/MCMC/SemanticEditors.hs new file mode 100644 --- /dev/null +++ b/MCMC/SemanticEditors.hs @@ -0,0 +1,74 @@ +-- | Semantic editor combinators +-- +-- + +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 diff --git a/MCMC/Tests.hs b/MCMC/Tests.hs new file mode 100644 --- /dev/null +++ b/MCMC/Tests.hs @@ -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)) + diff --git a/MCMC/Types.hs b/MCMC/Types.hs new file mode 100644 --- /dev/null +++ b/MCMC/Types.hs @@ -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) diff --git a/Tests.hs b/Tests.hs deleted file mode 100644 --- a/Tests.hs +++ /dev/null @@ -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] diff --git a/mcmc-samplers.cabal b/mcmc-samplers.cabal --- a/mcmc-samplers.cabal +++ b/mcmc-samplers.cabal @@ -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