diff --git a/Actions.hs b/Actions.hs new file mode 100644 --- /dev/null +++ b/Actions.hs @@ -0,0 +1,72 @@ +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 new file mode 100644 --- /dev/null +++ b/Distributions.hs @@ -0,0 +1,257 @@ +{-# 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/Gibbs.hs b/Gibbs.hs new file mode 100644 --- /dev/null +++ b/Gibbs.hs @@ -0,0 +1,195 @@ +-- | Gibbs sampling with a Naive Bayes model + +-- Details and notation based on: +-- http://www.cs.umd.edu/~hardisty/papers/gsfu.pdf + +import Control.Applicative +import Control.Monad +import Control.Monad.Primitive +import System.Random.MWC hiding (initialize) +import Statistics.Distribution +import Statistics.Distribution.Beta +import Statistics.Distribution.Gamma +import qualified Data.Vector.Unboxed as U +import qualified Data.Vector as V +import qualified Data.Map as M + +type HyperParams = U.Vector Int +type Theta = U.Vector Double -- Distribution over words in a document + +type Doc = [String] +type Label = Bool +type WordCounts = M.Map String Int + +data Point = Point { observe :: Doc, label :: Label } + deriving (Show) + +-- A point "augmented" with word counts information +data AugPoint = PnC { point :: Point, counts :: WordCounts } + deriving (Show) + +type Corpus = V.Vector AugPoint + +data Info = Info { dataSet :: Corpus + , thetas :: M.Map Label Theta + , wcInfo :: M.Map Label WordCounts + , numDocs :: M.Map Label Int } + +type Sample = V.Vector Label + +-- Beta(1,1), i.e, uniform +beta :: Gen RealWorld -> IO Double +beta = genContVar (betaDistr 1 1) + +-- Dirichlet as a vector of samples from Gamma +dirichlet :: HyperParams -> Gen RealWorld -> IO Theta +dirichlet hps gen = do + let gamma_draw hp = genContVar (gammaDistr (fromIntegral hp) 1) gen + ys <- U.mapM gamma_draw hps + let y_sum = U.sum ys + return $ U.map (/y_sum) ys + +bernoulli :: Double -> Gen RealWorld -> IO Bool +bernoulli p gen = uniform gen >>= return . (>p) + +tokens :: V.Vector String +tokens = V.fromList.words $ "I went to the square and saw a foolish \ + \politician who would not stop blabbering" + +capV :: Int +capV = V.length tokens + +capN :: Int +capN = 20 + +-- Docsize +kmax :: Int +kmax = 36 + +-- Generate a "bag of words" +gen_bag :: Gen RealWorld -> IO Doc +gen_bag g = let bagger b 0 = return b + bagger b n = do + i <- uniformR (0,capV-1) g + bagger ((tokens V.! i):b) $ n-1 + in uniformR (1,kmax) g >>= bagger [] + +doc_counts :: Point -> WordCounts +doc_counts (Point doc _) = foldr (M.adjust (+1)) zeroes doc + where ts = V.toList tokens + zeroes = M.fromList $ zip ts [0..] + +-- Assume an ordering of [a], in this case [true-value, false-value] +label_map :: [a] -> M.Map Label a +label_map vl = M.fromList $ zip [True, False] vl + +initialize :: Gen RealWorld -> IO Info +initialize g = do + p <- beta g + let gen_point = Point <$> gen_bag g <*> bernoulli p g + points <- V.fromList <$> replicateM capN gen_point + let corpus = V.zipWith PnC points $ V.map doc_counts points + (trues, falses) = V.unstablePartition (label.point) corpus + collect_counts ps = V.foldl1 (M.unionWith (+)) $ V.map counts ps + wcMap = label_map $ map collect_counts [trues, falses] + nums = label_map [V.length trues, V.length falses] + thets <- label_map <$> (replicateM 2 $ dirichlet (U.replicate capV 1) g) + return $ Info corpus thets wcMap nums + +cond_prob :: Int -> Theta -> WordCounts -> Double +cond_prob c_x theta_x wc_j = + let prod i t = (*) $ (^) t (wc_j M.! (tokens V.! i)) + p = U.ifoldr prod 1.0 theta_x + c = (/) (fromIntegral c_x) (fromIntegral $ capN + 1) + in c * p + +sample_label :: Int -> Info -> Gen RealWorld -> IO Label +sample_label j (Info dat thetas _ nums) gen = do + let wc_j = counts $ dat V.! j + pTrue = cond_prob (nums M.! True) (thetas M.! True) wc_j + pFalse = cond_prob (nums M.! False) (thetas M.! False) wc_j + pNorm = (/) pTrue $ pTrue + pFalse + bernoulli pNorm gen + +assign_label :: Int -> Label -> Info -> Info +assign_label j lab (Info d t w n) = + let ap = d V.! j + p = point ap + new_ap = PnC (Point (observe p) lab) (counts ap) + new_d = (V.//) d [(j,new_ap)] + in Info new_d t w n + +type WCUpdate = WordCounts -> WordCounts + +update_wc :: WCUpdate -> Label -> Info -> Info +update_wc fun lab (Info d t wc n) = + Info d t (M.adjust fun lab wc) n + +type NumUpdate = Int -> Int + +update_num :: NumUpdate -> Label -> Info -> Info +update_num fun lab (Info d t w nums) = + Info d t w (M.adjust fun lab nums) + +sampler :: Gen RealWorld -> Info -> IO Info +sampler gen info = do + let ind_ds = V.indexed $ dataSet info + f acc (j,ap) = do + let lab = (label . point) ap + sub_fun = M.differenceWith (\b a -> Just (a-b)) (counts ap) + pre_sample_info = update_num (flip (-) 1) lab $ + update_wc sub_fun lab acc + new_lab <- sample_label j pre_sample_info gen + let post_sample_info = assign_label j new_lab pre_sample_info + return $ update_wc (M.unionWith (+) (counts ap)) new_lab $ + update_num (+1) new_lab post_sample_info + V.foldM f info ind_ds -- TODO: Check whether foldM goes left-to-right + +new_thetas :: Gen RealWorld -> Info -> IO Info +new_thetas gen (Info d _ w n) = do + let f wc i = (+1) $ wc M.! (tokens V.! i) + hyperparams = M.map (U.generate capV) $ M.map f w + thetaT <- dirichlet (hyperparams M.! True) gen + thetaF <- dirichlet (hyperparams M.! False) gen + return $ Info d (label_map [thetaT, thetaF]) w n + +capT :: Int +capT = 1000 + +gibbs :: Gen RealWorld -> IO Sample +gibbs g = do + let loop 0 info = return info + loop t info = sampler g info >>= new_thetas g >>= loop (t-1) + info <- initialize g >>= loop capT + return $ V.map (label.point) (dataSet info) + +main :: IO () +main = testGibbs + +-- Tests + +testPRNG :: IO () +testPRNG = do + gen <- createSystemRandom + beta gen >>= print + beta gen >>= print + +testTheta :: IO () +testTheta = createSystemRandom >>= dirichlet (U.replicate capV 1) >>= print + +testBag :: IO () +testBag = createSystemRandom >>= gen_bag >>= print + +testInit :: IO () +testInit = createSystemRandom >>= initialize >>= print.(V.map counts).dataSet + +testCondProb :: IO () +testCondProb = do + info <- createSystemRandom >>= initialize + let pTrue = cond_prob (numDocs info M.! True) (thetas info M.! True) $ + counts $ dataSet info V.! 0 + putStrLn $ "P(L_0=True | Initial) = " ++ show pTrue + +testGibbs :: IO () +testGibbs = createSystemRandom >>= gibbs >>= print + diff --git a/Kernels.hs b/Kernels.hs new file mode 100644 --- /dev/null +++ b/Kernels.hs @@ -0,0 +1,103 @@ +{-# 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/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,24 @@ +This is free and unencumbered software released into the public domain. + +Anyone is free to copy, modify, publish, use, compile, sell, or +distribute this software, either in source code form or as a compiled +binary, for any purpose, commercial or non-commercial, and by any +means. + +In jurisdictions that recognize copyright laws, the author or authors +of this software dedicate any and all copyright interest in the +software to the public domain. We make this dedication for the benefit +of the public at large and to the detriment of our heirs and +successors. We intend this dedication to be an overt act of +relinquishment in perpetuity of all present and future rights to this +software under copyright law. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR +OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, +ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +OTHER DEALINGS IN THE SOFTWARE. + +For more information, please refer to diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,18 @@ +Samplers +======== + +### Here lies a library of combinators for MCMC kernels and proposals +- The relevant modules are `Kernels`, `Distributions`, and `Actions` +- See `Tests.hs` for some examples on how this library can be used +- Needs the [hmatrix](http://hackage.haskell.org/package/hmatrix) package + - Might need to do `cabal install hmatrix` + +##### On Gibbs.hs + +- The current implementation is for a Naive Bayes model +- TODO: + - Use an existing, "real" dataset instead of randomly generating sentences + - See which words appear most frequently for each label/class + - Average over all theta estimates and return top 10 and bottom 10 words + according to these averages + - Implement burn-in and lag (to decrease autocorrelation) diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/Tests.hs b/Tests.hs new file mode 100644 --- /dev/null +++ b/Tests.hs @@ -0,0 +1,88 @@ +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 new file mode 100644 --- /dev/null +++ b/mcmc-samplers.cabal @@ -0,0 +1,24 @@ +-- Initial mcmc-samplers.cabal generated by cabal init. For further +-- 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 +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 +cabal-version: >=1.10 + +library + exposed-modules: Actions, Distributions, Kernels + -- other-modules: + 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 + -- hs-source-dirs: + default-language: Haskell2010