mcmc-samplers (empty) → 0.1.0.0
raw patch · 9 files changed
+783/−0 lines, 9 filesdep +basedep +containersdep +hmatrixsetup-changed
Dependencies added: base, containers, hmatrix, mwc-random, primitive, statistics, vector
Files
- Actions.hs +72/−0
- Distributions.hs +257/−0
- Gibbs.hs +195/−0
- Kernels.hs +103/−0
- LICENSE +24/−0
- README.md +18/−0
- Setup.hs +2/−0
- Tests.hs +88/−0
- mcmc-samplers.cabal +24/−0
+ Actions.hs view
@@ -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+
+ Distributions.hs view
@@ -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
+ Gibbs.hs view
@@ -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+
+ Kernels.hs view
@@ -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+
+ LICENSE view
@@ -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 <http://unlicense.org/>
+ README.md view
@@ -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)
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ Tests.hs view
@@ -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]
+ mcmc-samplers.cabal view
@@ -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