We create a simple monad called PL which allows for a single operation: sampling
from a uniform distribution. We then exploit this to implement MCMC using metropolis hastings,
which is used to sample from arbitrary distributions. Bonus is a small library to render sparklines
in the CLI.
For next time:
- Using applicative to speed up computations by exploiting parallelism
- Conditioning of a distribution wrt a variable
§ Source code
{-# LANGUAGE GeneralizedNewtypeDeriving #-}{-# LANGUAGE GADTs #-}{-# LANGUAGE StandaloneDeriving #-}{-# LANGUAGE FlexibleContexts #-}{-# LANGUAGE FlexibleInstances #-}{-# LANGUAGE UndecidableInstances #-}{-# LANGUAGE DeriveFunctor #-}import System.Randomimport Data.List(sort, nub)import Data.Proxyimport Control.Monad (replicateM)import qualified Data.Map as M-- | Loop a monadic computation.mLoop :: Monad m => (a -> m a) -- ^ loop -> Int -- ^ number of times to run -> a -- initial value -> m a -- final valuemLoop _ 0 a = return amLoop f n a = f a >>= mLoop f (n - 1)-- | Utility library for drawing sparklines-- | List of characters that represent sparklinessparkchars :: Stringsparkchars = "_▁▂▃▄▅▆▇█"-- Convert an int to a sparkline characternum2spark :: RealFrac a => a -- ^ Max value -> a -- ^ Current value -> Charnum2spark maxv curv = sparkchars !! (floor $ (curv / maxv) * (fromIntegral (length sparkchars - 1)))series2spark :: RealFrac a => [a] -> Stringseries2spark vs = let maxv = if null vs then 0 else maximum vs in map (num2spark maxv) vsseriesPrintSpark :: RealFrac a => [a] -> IO ()seriesPrintSpark = putStrLn . series2spark-- Probabilities-- ============type F = Float-- | probability densitynewtype P = P { unP :: Float } deriving(Num)-- | prob. distributions over space anewtype D a = D { runD :: a -> P }uniform :: Int -> D auniform n = D $ \_ -> P $ 1.0 / (fromIntegral $ n)(>$<) :: Contravariant f => (b -> a) -> f a -> f b(>$<) = cofmapinstance Contravariant D where cofmap f (D d) = D (d . f)-- | Normal distribution with given meannormalD :: Float -> D FloatnormalD mu = D $ \f -> P $ exp (- ((f-mu)^2))-- | Distribution that takes on value x^p for 1 <= x <= 2. Is normalizedpolyD :: Float -> D FloatpolyD p = D $ \f -> P $ if 1 <= f && f <= 2 then (f ** p) * (p + 1) / (2 ** (p+1) - 1) else 0class Contravariant f where cofmap :: (b -> a) -> f a -> f bdata PL next where Ret :: next -> PL next -- ^ return a value Sample01 :: (Float -> PL next) -> PL next -- ^ sample uniformly from a [0, 1) distributioninstance Monad PL where return = Ret (Ret a) >>= f = f a (Sample01 float2plnext) >>= next2next' = Sample01 $ \f -> float2plnext f >>= next2next'instance Applicative PL where pure = return ff <*> fx = do f <- ff x <- fx return $ f xinstance Functor PL where fmap f plx = do x <- plx return $ f x-- | operation to sample from [0, 1)sample01 :: PL Floatsample01 = Sample01 Ret-- | Run one step of MH on a distribution to obtain a (correlated) samplemhStep :: (a -> Float) -- ^ function to score sample with, proportional to distribution -> (a -> PL a) -- ^ Proposal program -> a -- current sample -> PL amhStep f q a = do a' <- q a let alpha = f a' / f a -- acceptance ratio u <- sample01 return $ if u <= alpha then a' else a-- Typeclass that can provide me with data to run MCMC on itclass MCMC a where arbitrary :: a uniform2val :: Float -> ainstance MCMC Float where arbitrary = 0 -- map [0, 1) -> (-infty, infty) uniform2val v = tan (-pi/2 + pi * v){--- | Any enumerable object has a way to get me the starting point for MCMCinstance (Bounded a, Enum a) => MCMC a where arbitrary = toEnum 0 uniform2val v = let maxf = fromIntegral . fromEnum $ maxBound minf = fromIntegral . fromEnum $ minBound in toEnum $ floor $ minf + v * (maxf - minf)-}-- | Run MH to sample from a distributionmh :: (a -> Float) -- ^ function to score sample with -> (a -> PL a) -- ^ proposal program -> a -- ^ current sample -> PL amh f q a = mLoop (mhStep f q) 100 $ a-- | Construct a program to sample from an arbitrary distribution using MCMCmhD :: MCMC a => D a -> PL amhD (D d) = let scorer = (unP . d) proposal _ = do f <- sample01 return $ uniform2val f in mh scorer proposal arbitrary-- | Run the probabilistic value to get a samplesample :: RandomGen g => g -> PL a -> (a, g)sample g (Ret a) = (a, g)sample g (Sample01 f2plnext) = let (f, g') = random g in sample g' (f2plnext f)-- | Sample n values from the distributionsamples :: RandomGen g => Int -> g -> PL a -> ([a], g)samples 0 g _ = ([], g)samples n g pl = let (a, g') = sample g pl (as, g'') = samples (n - 1) g' pl in (a:as, g'')-- | count fraction of times value occurs in listoccurFrac :: (Eq a) => [a] -> a -> FloatoccurFrac as a = let noccur = length (filter (==a) as) n = length as in (fromIntegral noccur) / (fromIntegral n)-- | Produce a distribution from a PL by using the sampler to sample N timesdistribution :: (Eq a, Num a, RandomGen g) => Int -> g -> PL a -> (D a, g)distribution n g pl = let (as, g') = samples n g pl in (D (\a -> P (occurFrac as a)), g')-- | biased coincoin :: Float -> PL Int -- 1 with prob. p1, 0 with prob. (1 - p1)coin p1 = do Sample01 (\f -> Ret $ if f < p1 then 1 else 0)-- | Create a histogram from values.histogram :: Int -- ^ number of buckets -> [Float] -- values -> [Int]histogram nbuckets as = let minv :: Float minv = minimum as maxv :: Float maxv = maximum as -- value per bucket perbucket :: Float perbucket = (maxv - minv) / (fromIntegral nbuckets) bucket :: Float -> Int bucket v = floor (v / perbucket) bucketed :: M.Map Int Int bucketed = foldl (\m v -> M.insertWith (+) (bucket v) 1 m) mempty as in map snd . M.toList $ bucketedprintSamples :: (Real a, Eq a, Ord a, Show a) => String -> [a] -> IO ()printSamples s as = do putStrLn $ "***" <> s putStrLn $ " samples: " <> series2spark (map toRational as)printHistogram :: [Float] -> IO ()printHistogram samples = putStrLn $ series2spark (map fromIntegral . histogram 10 $ samples)-- | Given a coin bias, take samples and print biasprintCoin :: Float -> IO ()printCoin bias = do let g = mkStdGen 1 let (tosses, _) = samples 100 g (coin bias) printSamples ("bias: " <> show bias) tosses-- | Create normal distribution as sum of uniform distributions.normal :: PL Floatnormal = fromIntegral . sum <$> (replicateM 5 (coin 0.5))main :: IO ()main = do printCoin 0.01 printCoin 0.99 printCoin 0.5 printCoin 0.7 putStrLn $ "normal distribution using central limit theorem: " let g = mkStdGen 1 let (nsamples, _) = samples 1000 g normal -- printSamples "normal: " nsamples printHistogram nsamples putStrLn $ "normal distribution using MCMC: " let (mcmcsamples, _) = samples 1000 g (mhD $ normalD 0.5) printHistogram mcmcsamples putStrLn $ "sampling from x^4 with finite support" let (mcmcsamples, _) = samples 1000 g (mhD $ polyD 4) printHistogram mcmcsamples§ Output
***bias: 1.0e-2 samples: ________________________________________█_█________***bias: 0.99 samples: ███████████████████████████████████████████████████***bias: 0.5 samples: __█____█__███_███_█__█_█___█_█_██___████████__█_███***bias: 0.7 samples: __█__█_█__███_█████__███_█_█_█_██_█_████████__█████normal distribution using central limit theorem:_▄▇█▄_normal distribution using MCMC:__▁▄█▅▂▁___sampling from x^4 with finite support▁▁▃▃▃▄▅▆▇█_