diff --git a/Control/Monad/MC.hs b/Control/Monad/MC.hs deleted file mode 100644 --- a/Control/Monad/MC.hs +++ /dev/null @@ -1,14 +0,0 @@ ------------------------------------------------------------------------------ --- | --- Module : Control.Monad.MC --- Copyright : Copyright (c) , Patrick Perry --- License : BSD3 --- Maintainer : Patrick Perry --- Stability : experimental --- - -module Control.Monad.MC ( - module Control.Monad.MC.GSL - ) where - -import Control.Monad.MC.GSL diff --git a/Control/Monad/MC/GSL.hs b/Control/Monad/MC/GSL.hs deleted file mode 100644 --- a/Control/Monad/MC/GSL.hs +++ /dev/null @@ -1,271 +0,0 @@ -{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, UndecidableInstances #-} ------------------------------------------------------------------------------ --- | --- Module : Control.Monad.MC.GSL --- Copyright : Copyright (c) , Patrick Perry --- License : BSD3 --- Maintainer : Patrick Perry --- Stability : experimental --- - -module Control.Monad.MC.GSL ( - -- * The Monte Carlo monad - MC, - runMC, - evalMC, - execMC, - unsafeInterleaveMC, - - -- * The Monte Carlo monad transformer - MCT, - runMCT, - evalMCT, - execMCT, - liftMCT, - unsafeInterleaveMCT, - - -- | Pure random number generator creation - RNG, - mt19937, - - -- | Random distributions - uniform, - uniformInt, - normal, - poisson, - - ) where - -import Control.Monad ( liftM, MonadPlus(..) ) -import Control.Monad.Cont ( MonadCont(..) ) -import Control.Monad.Error ( MonadError(..) ) -import Control.Monad.Reader ( MonadReader(..) ) -import Control.Monad.State ( MonadState(..) ) -import Control.Monad.Writer ( MonadWriter(..) ) -import Control.Monad.Trans ( MonadTrans(..), MonadIO(..) ) -import Data.Word -import System.IO.Unsafe ( unsafePerformIO ) - -import GSL.Random.Gen hiding ( mt19937 ) -import qualified GSL.Random.Gen as Gen -import GSL.Random.Dist - --- | A Monte Carlo monad with an internal random number generator. -newtype MC a = MC (RNG -> (a,RNG)) - --- | Run this Monte Carlo monad with the given initial random number generator, --- getting the result and the new random number generator. -runMC :: MC a -> RNG -> (a, RNG) -runMC (MC g) r = - let r' = unsafePerformIO $ cloneRNG r - in r' `seq` g r' -{-# NOINLINE runMC #-} - --- | Evaluate this Monte Carlo monad and throw away the final random number --- generator. Very much like @fst@ composed with @runMC@. -evalMC :: MC a -> RNG -> a -evalMC g r = fst $ runMC g r - --- | Exicute this Monte Carlo monad and return the final random number --- generator. Very much like @snd@ composed with @runMC@. -execMC :: MC a -> RNG -> RNG -execMC g r = snd $ runMC g r - --- | Get the baton from the Monte Carlo monad without performing any --- computations. Useful but dangerous. -unsafeInterleaveMC :: MC a -> MC a -unsafeInterleaveMC (MC m) = MC $ \r -> let - (a,_) = m r - in (a,r) - - -instance Functor MC where - fmap f (MC m) = MC $ \r -> let - (a,r') = m r - in (f a, r') - -instance Monad MC where - return a = MC $ \r -> (a,r) - (MC m) >>= k = - MC $ \r -> let - (a, r') = m r - (MC m') = k a - in m' r' - -instance MonadState RNG MC where - get = MC $ getHelp - put r' = MC $ putHelp r' - -getHelp :: RNG -> (RNG,RNG) -getHelp r = unsafePerformIO $ do - r' <- cloneRNG r - r' `seq` return (r',r) -{-# NOINLINE getHelp #-} - -putHelp :: RNG -> RNG -> ((),RNG) -putHelp r' r = unsafePerformIO $ do - io <- copyRNG r r' - io `seq` return ((),r) -{-# NOINLINE putHelp #-} - --- | A parameterizable Monte Carlo monad for encapsulating an inner --- monad. -newtype MCT m a = MCT (RNG -> m (a,RNG)) - --- | Similar to 'runMC'. -runMCT :: (Monad m) => MCT m a -> RNG -> m (a,RNG) -runMCT (MCT g) r = - let r' = unsafePerformIO $ cloneRNG r - in r' `seq` g r' -{-# NOINLINE runMCT #-} - --- | Similar to 'evalMC'. -evalMCT :: (Monad m) => MCT m a -> RNG -> m a -evalMCT g r = do - ~(a,_) <- runMCT g r - return a - --- | Similar to 'execMC'. -execMCT :: (Monad m) => MCT m a -> RNG -> m RNG -execMCT g r = do - ~(_,r') <- runMCT g r - return r' - --- | Take a Monte Carlo computations and lift it to an MCT computation. -liftMCT :: (Monad m) => MC a -> MCT m a -liftMCT (MC m) = MCT $ return . m - --- | Similar to 'unsafeInterleaveMC'. -unsafeInterleaveMCT :: (Monad m) => MCT m a -> MCT m a -unsafeInterleaveMCT (MCT g) = MCT $ \r -> do - ~(a,_) <- g r - return (a,r) - -instance (Monad m) => Functor (MCT m) where - fmap f (MCT m) = MCT $ \r -> do - ~(x, r') <- m r - return (f x, r') - -instance (Monad m) => Monad (MCT m) where - return a = MCT $ \r -> return (a,r) - - (MCT m) >>= k = - MCT $ \r -> do - ~(a,r') <- m r - let (MCT m') = k a - m' r' - - fail str = MCT $ \_ -> fail str - -instance (MonadPlus m) => MonadPlus (MCT m) where - mzero = MCT $ \_ -> mzero - - (MCT m) `mplus` (MCT n) = - MCT $ \r -> - let r' = unsafePerformIO $ cloneRNG r - in r' `seq` (m r `mplus` n r') - -instance (Monad m) => MonadState RNG (MCT m) where - get = MCT $ return . getHelp - put r' = MCT $ return . (putHelp r') - -instance MonadTrans MCT where - lift m = MCT $ \r -> do - a <- m - return (a,r) - -instance (MonadCont m) => MonadCont (MCT m) where - callCC f = MCT $ \r -> - callCC $ \c -> - let (MCT m) = (f (\a -> MCT $ \r' -> c (a, r'))) - in m r - -instance (MonadError e m) => MonadError e (MCT m) where - throwError = lift . throwError - (MCT m) `catchError` h = MCT $ \r -> - m r `catchError` \e -> let (MCT m') = h e in m' r - -instance (MonadIO m) => MonadIO (MCT m) where - liftIO = lift . liftIO - -instance (MonadReader r m) => MonadReader r (MCT m) where - ask = lift ask - local f (MCT m) = MCT $ \r -> - local f (m r) - -instance (MonadState s m) => MonadState s (MCT m) where - get = lift get - put = lift . put - -instance (MonadWriter w m) => MonadWriter w (MCT m) where - tell = lift . tell - listen (MCT m) = MCT $ \r -> do - ~((a,r'),w) <- listen (m r) - return ((a,w),r') - pass (MCT m) = MCT $ \r -> pass $ do - ~((a,f),r') <- m r - return ((a,r'),f) - --- | Get a Mersenne Twister random number generator seeded with the given --- value. -mt19937 :: Word64 -> RNG -mt19937 s = unsafePerformIO $ do - r <- newRNG Gen.mt19937 - setSeed r s - return r -{-# NOINLINE mt19937 #-} - --- | @uniformInt n@ generates an integer uniformly in the range @[0,n-1]@. --- It is an error to call this function with a non-positive value. -uniformInt :: Int -> MC Int -uniformInt n = MC $ uniformIntHelp n - -uniformIntHelp :: Int -> RNG -> (Int,RNG) -uniformIntHelp n r = unsafePerformIO $ do - x <- getUniformInt r n - x `seq` return (x,r) - --- | @uniform a b@ generates a value uniformly distributed in @[a,b)@. -uniform :: Double -> Double -> MC Double -uniform a b = MC $ uniformHelp a b - -uniformHelp :: Double -> Double -> RNG -> (Double,RNG) -uniformHelp a b r = unsafePerformIO $ do - x <- getFlat r a b - x `seq` return (x,r) -{-# NOINLINE uniformHelp #-} - --- | @normal mu sigma@ generates a Normal random variable with mean --- @mu@ and standard deviation @sigma@. -normal :: Double -> Double -> MC Double -normal mu sigma = MC $ normalHelp mu sigma - -normalHelp :: Double -> Double -> RNG -> (Double,RNG) -normalHelp mu sigma r = unsafePerformIO $ do - x <- liftM (mu +) $ getGaussian r sigma - x `seq` return (x,r) -{-# NOINLINE normalHelp #-} - --- | @poisson mu@ generates a Poisson random variable with mean @mu@. -poisson :: Double -> MC Int -poisson mu = MC $ poissonHelp mu - -poissonHelp :: Double -> RNG -> (Int,RNG) -poissonHelp mu r = unsafePerformIO $ do - x <- getPoisson r mu - x `seq` return (x,r) -{-# NOINLINE poissonHelp #-} - - -{- - - -unifInt :: (Monad m) => Int -> MCT m Int -unifInt n = MCT $ unifInt' n - -unifInt' :: (Monad m) => Int -> RNG -> m (Int,RNG) -unifInt' n r = - unsafePerformIO $ do - i <- rngUnifInt r n - i `seq` (return . return) (i,r) --} diff --git a/examples/Pi.hs b/examples/Pi.hs --- a/examples/Pi.hs +++ b/examples/Pi.hs @@ -1,6 +1,7 @@ import Control.Monad.MC import Control.Monad +import Data.List( foldl' ) import System.Environment( getArgs ) import Text.Printf( printf ) @@ -13,27 +14,30 @@ inUnitCircle :: (Double,Double) -> Bool inUnitCircle (x,y) = x*x + y*y <= 1 --- | Generate @n@ points in the unit box and count how many are in the --- unit circle. -countInBox :: Int -> MC Int -countInBox n = do - xs <- replicateM n unitBox - return $ count inUnitCircle xs - --- | Count how many times the predicate is true -count :: (a -> Bool) -> [a] -> Int -count f = length . filter f +-- | Given a list of indicators, return the sample mean and standard +-- error. +average :: [Bool] -> (Double,Double) +average is = let + (t,n) = foldl' count (0,0) is + p = toDouble t / toDouble n + se = sqrt (p * (1 - p) / toDouble n) + in (p, se) + where + count (t,n) i = let + t' = if i then t+1 else t + n' = n+1 + in t' `seq` n' `seq` (t',n') + toDouble = realToFrac . toInteger + -- | Compute a Monte Carlo estimate of pi based on @n@ samples. Return -- the estimate and the standard error of the estimate. computePi :: Int -> MC (Double,Double) computePi n = do - m <- countInBox n - let p = toDouble m / toDouble n - se = sqrt (p * (1 - p) / toDouble n) - return (4*p, 4*se) - where - toDouble = realToFrac . toInteger + is <- liftM (map inUnitCircle) (unsafeInterleaveMC $ replicateM n unitBox) + let (mu ,se ) = average is + (mu',se') = (4*mu,4*se) + return (mu',se') -- | Given an estimate and standard error, produce a 99% confidence -- interval based on the Central Limit Theorem @@ -58,8 +62,9 @@ -- inverval coverage :: Int -> Int -> MC Int coverage r n = do - liftM (count id) $ replicateM r (covers n) - + liftM count $ replicateM r (covers n) + where + count = length . filter id main = do [n] <- map read `fmap` getArgs diff --git a/examples/Poker.hs b/examples/Poker.hs new file mode 100644 --- /dev/null +++ b/examples/Poker.hs @@ -0,0 +1,95 @@ +module Main where + +import Control.Monad +import Control.Monad.MC +import Data.List +import Data.Map( Map ) +import qualified Data.Map as Map +import System.Environment +import Text.Printf + +-- | Data types for representing cards. An Ace has 'number' equal to @1@. +-- Jack, Queen, and King have numbers @11@, @12@, and @13@, respectively. +data Suit = Club | Diamond | Heart | Spade deriving (Eq, Show) +data Card = Card { number :: Int + , suit :: Suit + } + deriving (Eq, Show) + +-- | The number values of the aces and face cards. +ace, jack, queen, king :: Int +ace = 1 +jack = 11 +queen = 12 +king = 13 + +-- | A type for the various poker hands. +data Hand = HighCard | Pair | TwoPair | ThreeOfAKind | Straight | Flush + | FullHouse | FourOfAKind | StraightFlush + deriving (Eq, Show, Ord) + +-- | Determine the hand corresponding to a list of five cards. +hand :: [Card] -> Hand +hand cs = + case matches of + [1,1,1,1,1] -> case undefined of + _ | isStraight && isFlush -> StraightFlush + _ | isFlush -> Flush + _ | isStraight -> Straight + _ | otherwise -> HighCard + [1,1,1,2] -> Pair + [1,2,2] -> TwoPair + [1,1,3] -> ThreeOfAKind + [2,3] -> FullHouse + [1,4] -> FourOfAKind + where + (x:xs) = (sort . map number) cs + (s:ss) = map suit cs + + isStraight | x == ace && xs == [ 10..king ] = True + | otherwise = xs == [ x+1..x+4 ] + + isFlush = all (== s) ss + + matches = (sort . map length . group) (x:xs) + + +-- | Get a list of cards that make up a 52-card deck. +deck :: [Card] +deck = [ Card i s | i <- [ 1..13 ], s <- [ Club, Diamond, Heart, Spade ] ] + +-- | Deal a five-card hand by choosing a random subset of the deck. +deal :: (MonadMC m) => m [Card] +deal = sampleSubset 5 52 deck + +-- | A type for storing the frequencies of the various hands. +type HandCounts = Map Hand Int + +-- | An empty frequency count. +emptyCounts :: HandCounts +emptyCounts = Map.empty + +-- | Update the count of the hand corresponding to a list of five cards. +updateCounts :: HandCounts -> [Card] -> HandCounts +updateCounts counts cs = Map.insertWith' (+) (hand cs) 1 counts + + +main = do + [reps] <- map read `fmap` getArgs + main' reps + +main' reps = + let seed = 0 + counts = repeatMCWith updateCounts emptyCounts reps deal + `evalMC` mt19937 seed in do + printf "\n" + printf " Hand Count Probability 99%% Interval \n" + printf "-------------------------------------------------------\n" + forM_ ((reverse . Map.toAscList) counts) $ \(h,c) -> + let n = fromIntegral reps :: Double + p = fromIntegral c / n + se = sqrt (p * (1 - p) / n) + delta = 2.575829 * se + (l,u) = (p-delta, p+delta) in + printf "%-13s %7d %.6f (%.6f,%.6f)\n" (show h) c p l u + printf "\n" diff --git a/examples/Sampling.hs b/examples/Sampling.hs new file mode 100644 --- /dev/null +++ b/examples/Sampling.hs @@ -0,0 +1,58 @@ +module Main + where + +import Control.Monad.MC +import Control.Monad +import Data.List( foldl' ) +import System.Environment( getArgs ) +import Text.Printf( printf ) + + +-- | Sample from a binomial distribution with the given parameters. +binomial :: (MonadMC m) => Int -> Double -> m Int +binomial n p = let + q = 1 - p + probs = map (\i -> (fromIntegral $ n `choose` i) * p^^i * q^^(n-i)) [0..n] + in sampleIntWithWeights probs (n+1) + +-- | Get a sample confidence interval for the mean after @reps@ replications of +-- a binomial with the given parameters. +binomialMean :: (MonadMC m) => Int -> Double -> Int -> m (Double,Double) +binomialMean n p reps = + liftM (sampleCI 0.95) $ repeatMC reps $ liftM fromIntegral (binomial n p) + +-- | Compute @reps@ 95% confidence intervals for the mean of an @(n,p)@ +-- binormal based on samples of the given size, and record the number +-- of intervals that contain the true mean. +coverage :: (MonadMC m) => Int -> Double -> Int -> Int -> m Int +coverage n p size reps = + repeatMCWith + (\c ci -> if mu `inInterval` ci then c+1 else c) + 0 + reps + (binomialMean n p size) + where + mu = fromIntegral n * p + x `inInterval` (l,h) = x > l && x < h + +main = do + [reps] <- map read `fmap` getArgs + main' reps + +main' reps = + let seed = 0 + n = 10 + p = 0.2 + size = 500 + c = evalMC (coverage n p size reps) $ mt19937 seed in + printf "\nOf %d 95%%-intervals, %d contain the true value.\n" reps c + + +--------------------------- Utility functions ----------------------------- + +factorial :: Int -> Int +factorial n | n <= 0 = 1 + | otherwise = n * factorial (n-1) + +choose :: Int -> Int -> Int +choose n k = factorial n `div` (factorial (n-k) * factorial k) diff --git a/lib/Control/Monad/MC.hs b/lib/Control/Monad/MC.hs new file mode 100644 --- /dev/null +++ b/lib/Control/Monad/MC.hs @@ -0,0 +1,14 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Control.Monad.MC +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- Stability : experimental +-- + +module Control.Monad.MC ( + module Control.Monad.MC.GSL + ) where + +import Control.Monad.MC.GSL diff --git a/lib/Control/Monad/MC/Base.hs b/lib/Control/Monad/MC/Base.hs new file mode 100644 --- /dev/null +++ b/lib/Control/Monad/MC/Base.hs @@ -0,0 +1,88 @@ +{-# LANGUAGE TypeFamilies, MultiParamTypeClasses, FlexibleContexts #-} +----------------------------------------------------------------------------- +-- | +-- Module : Control.Monad.MC.Base +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- Stability : experimental +-- + +module Control.Monad.MC.Base ( + -- * MonadMC type classes + MonadMC(..), + HasRNG(..), + + ) where + +import qualified Control.Monad.MC.GSLBase as GSL + +class HasRNG m where + -- | The random number generator type for the monad. + type RNG m + +class (Monad m, HasRNG m) => MonadMC m where + -- | Get the current random number generator. + getRNG :: m (RNG m) + + -- | Set the current random number generator. + setRNG :: RNG m -> m () + + -- | @uniform a b@ generates a value uniformly distributed in @[a,b)@. + uniform :: Double -> Double -> m Double + + -- | @uniformInt n@ generates an integer uniformly in the range @[0,n-1]@. + -- It is an error to call this function with a non-positive value. + uniformInt :: Int -> m Int + + -- | @normal mu sigma@ generates a Normal random variable with mean + -- @mu@ and standard deviation @sigma@. + normal :: Double -> Double -> m Double + + -- | @poisson mu@ generates a Poisson random variable with mean @mu@. + poisson :: Double -> m Int + + -- | Get the baton from the Monte Carlo monad without performing any + -- computations. Useful but dangerous. + unsafeInterleaveMC :: m a -> m a + + +------------------------------- Instances ----------------------------------- + +instance HasRNG GSL.MC where + type RNG GSL.MC = GSL.RNG + +instance MonadMC GSL.MC where + getRNG = GSL.getRNG + {-# INLINE getRNG #-} + setRNG = GSL.setRNG + {-# INLINE setRNG #-} + uniform = GSL.uniform + {-# INLINE uniform #-} + uniformInt = GSL.uniformInt + {-# INLINE uniformInt #-} + normal = GSL.normal + {-# INLINE normal #-} + poisson = GSL.poisson + {-# INLINE poisson #-} + unsafeInterleaveMC = GSL.unsafeInterleaveMC + {-# INLINE unsafeInterleaveMC #-} + +instance (Monad m) => HasRNG (GSL.MCT m) where + type RNG (GSL.MCT m) = GSL.RNG + +instance (Monad m) => MonadMC (GSL.MCT m) where + getRNG = GSL.liftMCT GSL.getRNG + {-# INLINE getRNG #-} + setRNG r = GSL.liftMCT $ GSL.setRNG r + {-# INLINE setRNG #-} + uniform a b = GSL.liftMCT $ GSL.uniform a b + {-# INLINE uniform #-} + uniformInt n = GSL.liftMCT $ GSL.uniformInt n + {-# INLINE uniformInt #-} + normal mu sigma = GSL.liftMCT $ GSL.normal mu sigma + {-# INLINE normal #-} + poisson mu = GSL.liftMCT $ GSL.poisson mu + {-# INLINE poisson #-} + unsafeInterleaveMC = GSL.unsafeInterleaveMCT + {-# INLINE unsafeInterleaveMC #-} diff --git a/lib/Control/Monad/MC/Class.hs b/lib/Control/Monad/MC/Class.hs new file mode 100644 --- /dev/null +++ b/lib/Control/Monad/MC/Class.hs @@ -0,0 +1,34 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Control.Monad.MC.Class +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- Stability : experimental +-- + +module Control.Monad.MC.Class ( + -- * The Monte Carlo monad type class + HasRNG(..), + MonadMC, + + -- * Getting and setting the generator + getRNG, + setRNG, + + -- * Random distributions + uniform, + uniformInt, + normal, + poisson, + + module Control.Monad.MC.Sample, + module Control.Monad.MC.Repeat, + + -- * Interleaving computations + unsafeInterleaveMC + ) where + +import Control.Monad.MC.Base +import Control.Monad.MC.Sample +import Control.Monad.MC.Repeat diff --git a/lib/Control/Monad/MC/GSL.hs b/lib/Control/Monad/MC/GSL.hs new file mode 100644 --- /dev/null +++ b/lib/Control/Monad/MC/GSL.hs @@ -0,0 +1,34 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Control.Monad.MC.GSL +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- Stability : experimental +-- + +module Control.Monad.MC.GSL ( + -- * The Monte Carlo monad + MC, + runMC, + evalMC, + execMC, + + -- * The Monte Carlo monad transformer + MCT, + runMCT, + evalMCT, + execMCT, + + -- * Pure random number generator creation + RNG, + mt19937, + + -- * Overloaded Monte Carlo monad interface + module Control.Monad.MC.Class, + + ) where + +import Control.Monad.MC.GSLBase ( MC, runMC, evalMC, execMC, + MCT, runMCT, evalMCT, execMCT, RNG, mt19937 ) +import Control.Monad.MC.Class hiding ( RNG ) diff --git a/lib/Control/Monad/MC/GSLBase.hs b/lib/Control/Monad/MC/GSLBase.hs new file mode 100644 --- /dev/null +++ b/lib/Control/Monad/MC/GSLBase.hs @@ -0,0 +1,285 @@ +{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, UndecidableInstances #-} +----------------------------------------------------------------------------- +-- | +-- Module : Control.Monad.MC.GSLBase +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- Stability : experimental +-- + +module Control.Monad.MC.GSLBase ( + -- * The Monte Carlo monad + MC(..), + runMC, + evalMC, + execMC, + unsafeInterleaveMC, + + -- * The Monte Carlo monad transformer + MCT(..), + runMCT, + evalMCT, + execMCT, + unsafeInterleaveMCT, + liftMCT, + + -- * Pure random number generator creation + RNG, + mt19937, + + -- * Getting and setting the random number generator + getRNG, + setRNG, + + -- * Random number distributions + uniform, + uniformInt, + normal, + poisson, + ) where + +import Control.Monad ( liftM, MonadPlus(..) ) +import Control.Monad.Cont ( MonadCont(..) ) +import Control.Monad.Error ( MonadError(..) ) +import Control.Monad.Reader ( MonadReader(..) ) +import Control.Monad.State ( MonadState(..) ) +import Control.Monad.Writer ( MonadWriter(..) ) +import Control.Monad.Trans ( MonadTrans(..), MonadIO(..) ) +import Data.Word +import System.IO.Unsafe ( unsafePerformIO ) + +import GSL.Random.Gen hiding ( mt19937 ) +import qualified GSL.Random.Gen as Gen +import GSL.Random.Dist + +-- | A Monte Carlo monad with an internal random number generator. +newtype MC a = MC (RNG -> (a,RNG)) + +-- | Run this Monte Carlo monad with the given initial random number generator, +-- getting the result and the new random number generator. +runMC :: MC a -> RNG -> (a, RNG) +runMC (MC g) r = + let r' = unsafePerformIO $ cloneRNG r + in r' `seq` g r' +{-# NOINLINE runMC #-} + +-- | Evaluate this Monte Carlo monad and throw away the final random number +-- generator. Very much like @fst@ composed with @runMC@. +evalMC :: MC a -> RNG -> a +evalMC g r = fst $ runMC g r + +-- | Exicute this Monte Carlo monad and return the final random number +-- generator. Very much like @snd@ composed with @runMC@. +execMC :: MC a -> RNG -> RNG +execMC g r = snd $ runMC g r + +unsafeInterleaveMC :: MC a -> MC a +unsafeInterleaveMC (MC m) = MC $ \r -> let + (a,_) = m r + in (a,r) + + +instance Functor MC where + fmap f (MC m) = MC $ \r -> let + (a,r') = m r + in (f a, r') + +instance Monad MC where + return a = MC $ \r -> (a,r) + {-# INLINE return #-} + + (MC m) >>= k = + MC $ \r -> let + (a, r') = m r + (MC m') = k a + in m' r' + {-# INLINE (>>=) #-} + +-- | A parameterizable Monte Carlo monad for encapsulating an inner +-- monad. +newtype MCT m a = MCT (RNG -> m (a,RNG)) + +-- | Similar to 'runMC'. +runMCT :: (Monad m) => MCT m a -> RNG -> m (a,RNG) +runMCT (MCT g) r = + let r' = unsafePerformIO $ cloneRNG r + in r' `seq` g r' +{-# NOINLINE runMCT #-} + +-- | Similar to 'evalMC'. +evalMCT :: (Monad m) => MCT m a -> RNG -> m a +evalMCT g r = do + ~(a,_) <- runMCT g r + return a + +-- | Similar to 'execMC'. +execMCT :: (Monad m) => MCT m a -> RNG -> m RNG +execMCT g r = do + ~(_,r') <- runMCT g r + return r' + +-- | Take a Monte Carlo computations and lift it to an MCT computation. +liftMCT :: (Monad m) => MC a -> MCT m a +liftMCT (MC m) = MCT $ return . m +{-# INLINE liftMCT #-} + +unsafeInterleaveMCT :: (Monad m) => MCT m a -> MCT m a +unsafeInterleaveMCT (MCT g) = MCT $ \r -> do + ~(a,_) <- g r + return (a,r) +{-# INLINE unsafeInterleaveMCT #-} + +instance (Monad m) => Functor (MCT m) where + fmap f (MCT m) = MCT $ \r -> do + ~(x, r') <- m r + return (f x, r') + {-# INLINE fmap #-} + +instance (Monad m) => Monad (MCT m) where + return a = MCT $ \r -> return (a,r) + {-# INLINE return #-} + + (MCT m) >>= k = + MCT $ \r -> do + ~(a,r') <- m r + let (MCT m') = k a + m' r' + {-# INLINE (>>=) #-} + + fail str = MCT $ \_ -> fail str + +instance (MonadPlus m) => MonadPlus (MCT m) where + mzero = MCT $ \_ -> mzero + {-# INLINE mzero #-} + + (MCT m) `mplus` (MCT n) = + MCT $ \r -> + let r' = unsafePerformIO $ cloneRNG r + in r' `seq` (m r `mplus` n r') + {-# NOINLINE mplus #-} + +instance MonadTrans MCT where + lift m = MCT $ \r -> do + a <- m + return (a,r) + {-# INLINE lift #-} + +instance (MonadCont m) => MonadCont (MCT m) where + callCC f = MCT $ \r -> + callCC $ \c -> + let (MCT m) = (f (\a -> MCT $ \r' -> c (a, r'))) + in m r + +instance (MonadError e m) => MonadError e (MCT m) where + throwError = lift . throwError + (MCT m) `catchError` h = MCT $ \r -> + m r `catchError` \e -> let (MCT m') = h e in m' r + +instance (MonadIO m) => MonadIO (MCT m) where + liftIO = lift . liftIO + {-# INLINE liftIO #-} + +instance (MonadReader r m) => MonadReader r (MCT m) where + ask = lift ask + local f (MCT m) = MCT $ \r -> + local f (m r) + +instance (MonadState s m) => MonadState s (MCT m) where + get = lift get + put = lift . put + +instance (MonadWriter w m) => MonadWriter w (MCT m) where + tell = lift . tell + listen (MCT m) = MCT $ \r -> do + ~((a,r'),w) <- listen (m r) + return ((a,w),r') + pass (MCT m) = MCT $ \r -> pass $ do + ~((a,f),r') <- m r + return ((a,r'),f) + +---------------------------- Random Number Generators ----------------------- + +getRNG :: MC RNG +getRNG = MC $ getHelp +{-# INLINE getRNG #-} + +getHelp :: RNG -> (RNG,RNG) +getHelp r = unsafePerformIO $ do + r' <- cloneRNG r + r' `seq` return (r',r) +{-# NOINLINE getHelp #-} + +setRNG :: RNG -> MC () +setRNG r' = MC $ setHelp r' +{-# INLINE setRNG #-} + +setHelp :: RNG -> RNG -> ((),RNG) +setHelp r' r = unsafePerformIO $ do + io <- copyRNG r r' + io `seq` return ((),r) +{-# NOINLINE setHelp #-} + +-- | Get a Mersenne Twister random number generator seeded with the given +-- value. +mt19937 :: Word64 -> RNG +mt19937 s = unsafePerformIO $ do + r <- newRNG Gen.mt19937 + setSeed r s + return r +{-# NOINLINE mt19937 #-} + + +-------------------------- Random Number Distributions ---------------------- + +uniform :: Double -> Double -> MC Double +uniform a b = MC $ uniformHelp a b +{-# INLINE uniform #-} + +uniformHelp :: Double -> Double -> RNG -> (Double,RNG) +uniformHelp 0 1 r = unsafePerformIO $ do + x <- getUniform r + x `seq` return (x,r) +uniformHelp a b r = unsafePerformIO $ do + x <- getFlat r a b + x `seq` return (x,r) +{-# NOINLINE uniformHelp #-} + +uniformInt :: Int -> MC Int +uniformInt n = MC $ uniformIntHelp n +{-# INLINE uniformInt #-} + +uniformIntHelp :: Int -> RNG -> (Int,RNG) +uniformIntHelp n r = unsafePerformIO $ do + x <- getUniformInt r n + x `seq` return (x,r) +{-# NOINLINE uniformIntHelp #-} + +normal :: Double -> Double -> MC Double +normal mu sigma = MC $ normalHelp mu sigma +{-# INLINE normal #-} + +normalHelp :: Double -> Double -> RNG -> (Double,RNG) +normalHelp 0 1 r = unsafePerformIO $ do + x <- getUGaussianRatioMethod r + x `seq` return (x,r) +normalHelp mu 1 r = unsafePerformIO $ do + x <- liftM (mu +) $ getUGaussianRatioMethod r + x `seq` return (x,r) +normalHelp 0 sigma r = unsafePerformIO $ do + x <- getGaussianRatioMethod r sigma + x `seq` return (x,r) +normalHelp mu sigma r = unsafePerformIO $ do + x <- liftM (mu +) $ getGaussianRatioMethod r sigma + x `seq` return (x,r) +{-# NOINLINE normalHelp #-} + +poisson :: Double -> MC Int +poisson mu = MC $ poissonHelp mu +{-# INLINE poisson #-} + +poissonHelp :: Double -> RNG -> (Int,RNG) +poissonHelp mu r = unsafePerformIO $ do + x <- getPoisson r mu + x `seq` return (x,r) +{-# NOINLINE poissonHelp #-} diff --git a/lib/Control/Monad/MC/Repeat.hs b/lib/Control/Monad/MC/Repeat.hs new file mode 100644 --- /dev/null +++ b/lib/Control/Monad/MC/Repeat.hs @@ -0,0 +1,54 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Control.Monad.MC.Repeat +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- Stability : experimental +-- + +module Control.Monad.MC.Repeat ( + -- * Averaging functions + repeatMC, + repeatMCWith, + + module Control.Monad.MC.Summary, + ) where + +import Control.Monad +import Control.Monad.MC.Base +import Control.Monad.MC.Summary +import Data.List( foldl' ) + +-- | Repeat a Monte Carlo generator the given number of times and return +-- the sample summary statistics. Note that this only works with +-- @Double@s. +repeatMC :: (MonadMC m) + => Int + -> m Double + -> m Summary +repeatMC = repeatMCWith update summary +{-# INLINE repeatMC #-} + +-- | Generalized version of 'repeatMC'. Run a Monte Carlo generator +-- the given number of times and accumulate the results. The accumulator +-- is strictly evaluated. +repeatMCWith :: (MonadMC m) + => (a -> b -> a) -- ^ accumulator + -> a -- ^ initial value + -> Int -- ^ number of repetitions + -> m b -- ^ generator + -> m a +repeatMCWith f a n mb = do + bs <- interleaveSequence $ replicate n mb + return $! foldl' f a bs +{-# INLINE repeatMCWith #-} + + +interleaveSequence :: (MonadMC m) => [m a] -> m [a] +interleaveSequence [] = return [] +interleaveSequence (m:ms) = unsafeInterleaveMC $ do + a <- m + as <- interleaveSequence ms + return (a:as) +{-# INLINE interleaveSequence #-} diff --git a/lib/Control/Monad/MC/Sample.hs b/lib/Control/Monad/MC/Sample.hs new file mode 100644 --- /dev/null +++ b/lib/Control/Monad/MC/Sample.hs @@ -0,0 +1,189 @@ +{-# LANGUAGE ScopedTypeVariables #-} +----------------------------------------------------------------------------- +-- | +-- Module : Control.Monad.MC.Sample +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- Stability : experimental +-- + +module Control.Monad.MC.Sample ( + -- * Sampling from lists + sample, + sampleWithWeights, + sampleSubset, + + -- * Sampling @Int@s + sampleInt, + sampleIntWithWeights, + sampleIntSubset, + + -- * Shuffling + shuffle, + shuffleInt, + ) where + +import Control.Monad +import Control.Monad.ST +import Control.Monad.MC.Base +import Control.Monad.MC.Walker + +import Data.Array.Base +import Data.Array.IArray +import Data.Array.ST +import Data.Array.Vector + +-- | @sample n xs@ samples a value uniformly from @take n xs@. The results +-- are undefined if @length xs@ is less than @n@. +sample :: (MonadMC m) => Int -> [a] -> m a +sample n xs = + sampleHelp n xs $ sampleInt n +{-# INLINE sample #-} + +-- | @sampleWithWeights ws n xs@ samples a value from @take n xs@, putting +-- weight @ws !! i@ on element @xs !! i@. The results +-- are undefined if @length xs@ or @length ws@ is less than @n@. +sampleWithWeights :: (MonadMC m) => [Double] -> Int -> [a] -> m a +sampleWithWeights ws n xs = + sampleHelp n xs $ sampleIntWithWeights ws n +{-# INLINE sampleWithWeights #-} + +-- | @sampleSubset k n xs@ samples a subset of size @k@ from @take n xs@ by +-- sampling without replacement. The return value is a list of length @k@ +-- with the elements in the subset in the order that they were sampled. Note +-- also that the elements are lazily generated. The results are undefined +-- if @k > n@ or if @length xs < n@. +sampleSubset :: (MonadMC m) => Int -> Int -> [a] -> m [a] +sampleSubset k n xs = + sampleListHelp n xs $ sampleIntSubset k n +{-# INLINE sampleSubset #-} + +sampleHelp :: (Monad m) => Int -> [a] -> m Int -> m a +sampleHelp n (xs :: [a]) f = let + arr = listArray (0,n-1) xs :: Array Int a + in liftM (unsafeAt arr) f + +sampleHelpUA :: (UA a, Monad m) => Int -> [a] -> m Int -> m a +sampleHelpUA n xs f = let + arr = newU n (\marr -> zipWithM_ (writeMU marr) [0..n-1] xs) + in liftM (indexU arr) f + +{-# RULES "sampleHelp/Double" forall n xs f. + sampleHelp n (xs :: [Double]) f = sampleHelpUA n xs f #-} +{-# RULES "sampleHelp/Int" forall n xs f. + sampleHelp n (xs :: [Int]) f = sampleHelpUA n xs f #-} + +sampleListHelp :: (Monad m) => Int -> [a] -> m [Int] -> m [a] +sampleListHelp n (xs :: [a]) f = let + arr = listArray (0,n-1) xs :: Array Int a + in liftM (map $ unsafeAt arr) f + +sampleListHelpUA :: (UA a, Monad m) => Int -> [a] -> m [Int] -> m [a] +sampleListHelpUA n xs f = let + arr = newU n (\marr -> zipWithM_ (writeMU marr) [0..n-1] xs) + in liftM (map $ indexU arr) f + +{-# RULES "sampleListHelp/Double" forall n xs f. + sampleListHelp n (xs :: [Double]) f = sampleListHelpUA n xs f #-} +{-# RULES "sampleListHelp/Int" forall n xs f. + sampleListHelp n (xs :: [Int]) f = sampleListHelpUA n xs f #-} + +-- | @sampleInt n@ samples integers uniformly from @[ 0..n-1 ]@. It is an +-- error to call this function with a non-positive @n@. +sampleInt :: (MonadMC m) => Int -> m Int +sampleInt n | n < 1 = fail "invalid argument" + | otherwise = uniformInt n +{-# INLINE sampleInt #-} + +-- | @sampleIntWithWeights ws n@ samples integers from @[ 0..n-1 ]@ with the +-- probability of choosing @i@ proportional to @ws !! i@. The list @ws@ must +-- have length equal to @n@. Also, the elements of @ws@ must be non-negative +-- with at least one nonzero entry. +sampleIntWithWeights :: (MonadMC m) => [Double] -> Int -> m Int +sampleIntWithWeights ws n = + let qjs = computeTable n ws + in liftM (indexTable qjs) (uniform 0 1) +{-# INLINE sampleIntWithWeights #-} + +-- | @sampleIntSubset k n@ samples a subset of size @k@ by sampling without +-- replacement from the integers @{ 0, ..., n-1 }@. The return value is a +-- list of length @k@ with the elements in the subset in the order that they +-- were sampled. Note also that the elements are lazily generated. +sampleIntSubset :: (MonadMC m) => Int -> Int -> m [Int] +sampleIntSubset k n | k < 0 = fail "negative subset size" + | k > n = fail "subset size is too big" + | otherwise = do + us <- randomIndices k n + return $ runST $ do + ints <- newMU n + sequence_ [ writeMU ints i i | i <- [0 .. n-1] ] + sampleIntSubsetHelp ints us (n-1) + where + randomIndices k' n' | k' == 0 = return [] + | otherwise = unsafeInterleaveMC $ do + u <- uniformInt n' + us <- randomIndices (k'-1) (n'-1) + return (u:us) + + sampleIntSubsetHelp _ [] _ = return [] + sampleIntSubsetHelp ints (u:us) n' = unsafeInterleaveST $ do + i <- readMU ints u + writeMU ints u =<< readMU ints n' + is <- sampleIntSubsetHelp ints us (n'-1) + return (i:is) +{-# INLINE sampleIntSubset #-} + +-- | @shuffle n xs@ randomly permutes the list @take n xs@ and returns +-- the result. All permutations of the elements of @xs@ are equally +-- likely. The results are undefined if @length xs@ is less than @n@. +shuffle :: (MonadMC m) => Int -> [a] -> m [a] +shuffle n (xs :: [a]) = + shuffleInt n >>= \swaps -> (return . runST) $ do + marr <- newListArray (0,n-1) xs :: ST s (STArray s Int a) + mapM_ (swap marr) swaps + getElems marr + where + swap marr (i,j) | i == j = return () + | otherwise = do + x <- unsafeRead marr i + y <- unsafeRead marr j + unsafeWrite marr i y + unsafeWrite marr j x +{-# INLINE shuffle #-} + +shuffleUA :: (UA a, MonadMC m) => Int -> [a] -> m [a] +shuffleUA n (xs :: [a]) = + shuffleInt n >>= \swaps -> (return . runST) $ do + marr <- newMU n + zipWithM_ (writeMU marr) [0 .. n-1] xs + mapM_ (swap marr) swaps + arr <- unsafeFreezeAllMU marr + return $ fromU arr + where + swap marr (i,j) | i == j = return () + | otherwise = do + x <- readMU marr i + y <- readMU marr j + writeMU marr i y + writeMU marr j x +{-# INLINE shuffleUA #-} + +{-# RULES "shuffle/Double" forall n xs. + shuffle n (xs :: [Double]) = shuffleUA n xs #-} +{-# RULES "shuffle/Int" forall n xs. + shuffle n (xs :: [Int]) = shuffleUA n xs #-} + + +-- | @shuffleInt n@ generates a sequence of swaps equivalent to a +-- uniformly-chosen random permutatation of the integers @{0, ..., n-1}@. +-- For an input of @n@, there are @n-1@ swaps, which are lazily generated. +shuffleInt :: (MonadMC m) => Int -> m [(Int,Int)] +shuffleInt n = + let shuffleIntHelp i | i <= 1 = return [] + | otherwise = unsafeInterleaveMC $ do + j <- uniformInt i + ijs <- shuffleIntHelp (i-1) + return $ (i-1,j):ijs in + shuffleIntHelp n +{-# INLINE shuffleInt #-} diff --git a/lib/Control/Monad/MC/Summary.hs b/lib/Control/Monad/MC/Summary.hs new file mode 100644 --- /dev/null +++ b/lib/Control/Monad/MC/Summary.hs @@ -0,0 +1,96 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Control.Monad.MC.Summary +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- Stability : experimental +-- + +module Control.Monad.MC.Summary ( + -- * Summary statistics + -- ** The @Summary@ data type + Summary, + summary, + update, + + -- ** @Summary@ properties + sampleSize, + sampleMean, + sampleVar, + sampleSD, + sampleSE, + sampleCI, + sampleMin, + sampleMax, + + ) where + +import GSL.Random.Dist( ugaussianPInv ) + +-- | A type for storing summary statistics for a data set including +-- sample size, min and max values, and first and second moments. +data Summary = S {-# UNPACK #-} !Int -- sample size + {-# UNPACK #-} !Double -- sample mean + {-# UNPACK #-} !Double -- sum of squares + {-# UNPACK #-} !Double -- sample min + {-# UNPACK #-} !Double -- sample max + +-- | Get an empty summary. +summary :: Summary +summary = S 0 0 0 (1/0) (-1/0) + +-- | Update the summary with a data point. +-- Running mean and variance computed as in Knuth, Vol 2, page 232, +-- 3rd edition, see http://www.johndcook.com/standard_deviation.html for +-- a description. +update :: Summary -> Double -> Summary +update (S n m s l h) x = + let n' = n+1 + delta = x - m + m' = m + delta / fromIntegral n' + s' = s + delta*(x - m') + l' = if x < l then x else l + h' = if x > h then x else h + in S n' m' s' l' h' + +-- | Get the sample size. +sampleSize :: Summary -> Int +sampleSize (S n _ _ _ _) = n + +-- | Get the sample mean. +sampleMean :: Summary -> Double +sampleMean (S _ m _ _ _) = m + +-- | Get the sample variance. +sampleVar :: Summary -> Double +sampleVar (S n _ s _ _) = s / fromIntegral (n - 1) + +-- | Get the sample standard deviation. +sampleSD :: Summary -> Double +sampleSD s = sqrt (sampleVar s) + +-- | Get the sample standard error. +sampleSE :: Summary -> Double +sampleSE s = sqrt (sampleVar s / fromIntegral (sampleSize s)) + +-- | Get a Central Limit Theorem-based confidence interval for the mean +-- with the specified coverage level. The level must be in the range @(0,1)@. +sampleCI :: Double -> Summary -> (Double,Double) +sampleCI level s | not (level > 0 && level < 1) = + error "level must be between 0 and 1" + | otherwise = + let alpha = (0.5 - level) + 0.5 + z = -(ugaussianPInv (0.5*alpha)) + se = sampleSE s + delta = z*se + xbar = sampleMean s + in (xbar-delta, xbar+delta) + +-- | Get the minimum of the sample. +sampleMin :: Summary -> Double +sampleMin (S _ _ _ l _) = l + +-- | Get the maximum of the sample. +sampleMax :: Summary -> Double +sampleMax (S _ _ _ _ h) = h diff --git a/lib/Control/Monad/MC/Walker.hs b/lib/Control/Monad/MC/Walker.hs new file mode 100644 --- /dev/null +++ b/lib/Control/Monad/MC/Walker.hs @@ -0,0 +1,159 @@ +{-# LANGUAGE TypeOperators #-} +----------------------------------------------------------------------------- +-- | +-- Module : Control.Monad.MC.Walker +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- Stability : experimental +-- +-- An implementation of Walker's Alias method for sampling from discrete +-- distributions. See section III.4 of Luc Devroye's book +-- "Non-Uniform Random Variate Generation", which is available on his +-- homepage, for a description of how it works. +module Control.Monad.MC.Walker ( + Table, + computeTable, + indexTable, + tableSize, + component, + ) where + +import Control.Monad +import Control.Monad.ST +import Data.Array.Vector + +-- | The table, which represents an equiprobable mixture of two-point +-- distributions. The @l@th entry of the table represents a mixture +-- distribution with weight @q[l]@ on @l@ and weight @(1-q[l])@ on @j[l]@. +-- The @l@th element of the table stores the pair @q[l] :*: j[l]@. +newtype Table = T (UArr (Double :*: Int)) + +-- | Get the @i@th mixture component. That is, return @q[i]@ and @j[i]@, +-- where the @i@th mixture component puts mass @q[i]@ on @i@ and mass +-- @1 - q[i]@ on @j[i]@. +component :: Table -> Int -> (Double,Int) +component (T qjs) i = let + (q' :*: j) = indexU qjs i + q = q' - fromIntegral i + in (q,j) + +-- | Compute the table for use in Walker's aliasing method. +computeTable :: Int -> [Double] -> Table +computeTable n ws = runST $ do + (qjs, sets) <- initTable n ws + breakLarger qjs sets + scaleTable qjs + liftM T $ unsafeFreezeAllMU qjs + +-- | Given an alias table and a number in the range [0,1), +-- get the corresponding sample in the table. +indexTable :: Table -> Double -> Int +indexTable (T qjs) u = let + n = lengthU qjs + nu = u * fromIntegral n + l = floor nu + (ql :*: jl) = indexU qjs l + in if nu < ql then l else jl + +-- | Get the size of the table +tableSize :: Table -> Int +tableSize (T qjs) = lengthU qjs + +-- | An intermediate result for use in computing a Table. +type STTable s = MUArr (Double :*: Int) s + +-- | A partition of indices into the sets /Greater/ and /Smaller/. The +-- indices of the /Smaller/ set are stored in positions @0, ..., numSmall - 1@, +-- and the indices of the /Greater/ set are stored in positions +-- @numSmall, ..., n-1@, where @n@ is the size of the underlying array. +data STPartition s = P !(MUArr Int s) + !Int + +-- | Given a list of weights, @ws@, compute corresponding probabilities, @ps@, +-- and store @map (n*) ps@ in the @qs@ array. Partition the probabilities +-- into two sets, /Greater/, and /Smaller/ based on whether or not +-- @q >= 1@ or @q < 1@. +initTable :: Int -> [Double] -> ST s (STTable s, STPartition s) +initTable n ws = do + when (n < 0) $ fail "negative table size" + sets <- newMU n :: ST s (MUArr Int s) + qjs <- newMU n :: ST s (MUArr (Double :*: Int) s) + + -- Store the weights in the table and compute their total. + total <- + foldM (\current (i,w) -> do + if w >= 0 + then do + writeMU qjs i (w :*: i) + return $! current + w + else + fail $ "negative probability" ) + 0 + (zip [0 .. n-1] ws) + + when (total == 0) $ fail "no positive probabilities given" + + -- scale the weights to get the qs, and partition the probabilites + -- into the two sets + let scale = fromIntegral n / total + nsmall <- liftM fst $ + foldM (\(smaller,greater) i -> do + p <- liftM fstS $ readMU qjs i + let q = scale*p + writeMU qjs i (q :*: i) + if q < 1 + then do + writeMU sets smaller i + return (smaller+1,greater) + else do + writeMU sets greater i + return (smaller,greater-1) ) + (0,n-1) + [0 .. n-1] + + return $ (qjs, P sets nsmall) + + +-- Given an initialized table and partition, compute the two-point +-- distributions by splitting the larger probabilites sccross multiple +-- distribions. +breakLarger :: STTable s -> STPartition s -> ST s () +breakLarger qjs (P sets nsmall) | nsmall == 0 = return () + | otherwise = let + n = lengthMU qjs + breakLargerHelp nsmall' i | nsmall' == n = return () + | i == n = return () + | otherwise = do + -- while Greater is not empty + -- choose k from Greater, l from Smaller + k <- readMU sets $ nsmall' + l <- readMU sets $ i + qk <- liftM fstS $ readMU qjs k + ql <- liftM fstS $ readMU qjs l + + -- set jl := k, finalize (ql,jl) + let jl = k + writeMU qjs l (ql :*: jl) + + -- set qk := qk - (1-ql) + let qk' = qk - (1-ql) + writeMU qjs k (qk' :*: k) + + -- if qk' < 1, move k from Greater to Smaller + let nsmall'' = if qk' < 1 then nsmall'+1 else nsmall' + + breakLargerHelp nsmall'' (i+1) + in + breakLargerHelp nsmall 0 + +-- Scale the probabilities in the table so that the lth entry +-- stores q[l] + l instead of q[l]. This helps when we are sampling +-- from the table. +scaleTable :: STTable s -> ST s () +scaleTable qjs = let + n = lengthMU qjs in + forM_ [ 0..(n-1) ] $ \l -> do + (ql :*: jl) <- readMU qjs l + writeMU qjs l ((ql + fromIntegral l) :*: jl) + diff --git a/monte-carlo.cabal b/monte-carlo.cabal --- a/monte-carlo.cabal +++ b/monte-carlo.cabal @@ -1,32 +1,48 @@ -name: monte-carlo -version: 0.1 -homepage: http://stat.stanford.edu/~patperry/code/monte-carlo -synopsis: A monad and transformer for Monte Carlo calculations. -description: - A monad and transformer for Monte Carlo calculations. The monads - carry and provide access to a random number generator. Importantly, - they only keep one copy of the generator state, and so are much more - efficient than MonadRandom. Currently, only the generator that comes - with the GNU Scientific Library is supported. - . -category: Math -license: BSD3 -license-file: LICENSE -copyright: (c) 2008. Patrick Perry -author: Patrick Perry -maintainer: Patrick Perry -cabal-version: >= 1.2.0 -build-type: Simple -tested-with: GHC ==6.8.3 - -extra-source-files: examples/Pi.hs +name: monte-carlo +version: 0.2 +license: BSD3 +license-file: LICENSE +author: Patrick Perry +maintainer: Patrick Perry +homepage: http://github.com/patperry/monte-carlo +category: Math +synopsis: A monad and transformer for Monte Carlo calculations. +description: A monad and transformer for Monte Carlo calculations. The + monads carry and provide access to a random number generator. + Importantly, they only keep one copy of the generator state, + and so are much more efficient than MonadRandom. Currently, + only the generator that comes with the GNU Scientific Library + (GSL) is supported. +build-type: Simple +stability: experimental +cabal-version: >= 1.2.3 +extra-source-files: examples/Pi.hs, examples/Sampling.hs examples/Poker.hs + tests/Main.hs tests/Makefile library - exposed-modules: Control.Monad.MC - Control.Monad.MC.GSL - - ghc-options: -Wall - extensions: MultiParamTypeClasses, FlexibleInstances, - UndecidableInstances + build-depends: array, base, mtl, gsl-random >=0.2.3, uvector + + exposed-modules: + Control.Monad.MC + Control.Monad.MC.Class + + other-modules: + Control.Monad.MC.Base + Control.Monad.MC.GSL + Control.Monad.MC.GSLBase + Control.Monad.MC.Repeat + Control.Monad.MC.Sample + Control.Monad.MC.Summary + Control.Monad.MC.Walker + + extensions: + FlexibleContexts, + FlexibleInstances, + MultiParamTypeClasses, + ScopedTypeVariables, + TypeFamilies, + TypeOperators, + UndecidableInstances - build-depends: base, mtl, gsl-random + hs-source-dirs: lib + ghc-options: -Wall diff --git a/tests/Main.hs b/tests/Main.hs new file mode 100644 --- /dev/null +++ b/tests/Main.hs @@ -0,0 +1,172 @@ + +module Main where + +import Debug.Trace +import Control.Monad +import Data.AEq +import Data.List +import System.IO +import System.Environment +import System.Random +import Text.Printf +import Test.QuickCheck + +import Control.Monad.MC.Walker + + +prop_table_probs (Weights n ws) = + let table = computeTable n ws + in all (\i -> probOf table i ~== ps !! i) [0..n-1] + where + ps = probsFromWeights ws + +prop_table_index (Weights n ws) (Unif u) = + let table = computeTable n ws + i = indexTable table u + in i >= 0 && i < n && (ws !! i > 0) + +tests_Walker = [ ("table probabilities", mytest prop_table_probs) + , ("table indexing" , mytest prop_table_index) + ] + +probOf table i = + (((sum . map ((1-) . fst) . filter ((==i) . snd)) + (map (component table) [0..n-1])) + + (fst . component table) i) / fromIntegral n + where + n = tableSize table + +------------------------------- Utility functions --------------------------- + +probsFromWeights ws = let + w = sum ws + ps = map (/w) ws + in ps + +------------------------------- Test generators ----------------------------- + +posInt :: Gen Int +posInt = do + n <- arbitrary + return $! abs n + 1 + +weight :: Gen Double +weight = do + w <- liftM abs arbitrary + if w < infty then return w else weight + where + infty = 1/0 + +weights :: Int -> Gen [Double] +weights n = do + ws <- replicateM n weight + if not (all (== 0) ws) then return ws else return $ replicate n 1.0 + +unif :: Gen Double +unif = do + u <- choose (0,1) + if u == 1 then return 0 else return u + +data Weights = Weights Int [Double] deriving Show +instance Arbitrary Weights where + arbitrary = do + n <- posInt + ws <- weights n + return $ Weights n ws + + coarbitrary (Weights n ws) = + coarbitrary (n,ws) + +data Unif = Unif Double deriving Show +instance Arbitrary Unif where + arbitrary = liftM Unif unif + coarbitrary (Unif u) = coarbitrary u + +------------------------------------------------------------------------ +-- +-- QC driver ( taken from xmonad-0.6 ) +-- + +debug = False + +mytest :: Testable a => a -> Int -> IO (Bool, Int) +mytest a n = mycheck defaultConfig + { configMaxTest=n + , configEvery = \n args -> let s = show n in s ++ [ '\b' | _ <- s ] } a + -- , configEvery= \n args -> if debug then show n ++ ":\n" ++ unlines args else [] } a + +mycheck :: Testable a => Config -> a -> IO (Bool, Int) +mycheck config a = do + rnd <- newStdGen + mytests config (evaluate a) rnd 0 0 [] + +mytests :: Config -> Gen Result -> StdGen -> Int -> Int -> [[String]] -> IO (Bool, Int) +mytests config gen rnd0 ntest nfail stamps + | ntest == configMaxTest config = done "OK," ntest stamps >> return (True, ntest) + | nfail == configMaxFail config = done "Arguments exhausted after" ntest stamps >> return (True, ntest) + | otherwise = + do putStr (configEvery config ntest (arguments result)) >> hFlush stdout + case ok result of + Nothing -> + mytests config gen rnd1 ntest (nfail+1) stamps + Just True -> + mytests config gen rnd1 (ntest+1) nfail (stamp result:stamps) + Just False -> + putStr ( "Falsifiable after " + ++ show ntest + ++ " tests:\n" + ++ unlines (arguments result) + ) >> hFlush stdout >> return (False, ntest) + where + result = generate (configSize config ntest) rnd2 gen + (rnd1,rnd2) = split rnd0 + +done :: String -> Int -> [[String]] -> IO () +done mesg ntest stamps = putStr ( mesg ++ " " ++ show ntest ++ " tests" ++ table ) + where + table = display + . map entry + . reverse + . sort + . map pairLength + . group + . sort + . filter (not . null) + $ stamps + + display [] = ".\n" + display [x] = " (" ++ x ++ ").\n" + display xs = ".\n" ++ unlines (map (++ ".") xs) + + pairLength xss@(xs:_) = (length xss, xs) + entry (n, xs) = percentage n ntest + ++ " " + ++ concat (intersperse ", " xs) + + percentage n m = show ((100 * n) `div` m) ++ "%" + +------------------------------------------------------------------------ + + + +main :: IO () +main = do + args <- getArgs + let n = if null args then 100 else read (head args) + + (results, passed) <- liftM unzip $ + foldM ( \prev (name,subtests) -> do + printf "\n%s\n" name + printf "%s\n" $ replicate (length name) '-' + cur <- mapM (\(s,a) -> printf "%-30s: " s >> a n) subtests + return (prev ++ cur) + ) + [] + tests + + printf "\nPassed %d tests!\n\n" (sum passed) + when (not . and $ results) $ fail "\nNot all tests passed!" + where + + tests = [ ("Walker" , tests_Walker) + ] diff --git a/tests/Makefile b/tests/Makefile new file mode 100644 --- /dev/null +++ b/tests/Makefile @@ -0,0 +1,17 @@ +all: + ghc -O -i. -i../lib Main.hs --make -o test-mc + ./test-mc + +hpc: + ghc -fforce-recomp -i. -i../lib -fhpc --make Main.hs -o test-mc + rm -f test-mc.tix + ./test-mc + hpc markup test-mc + +clean: + find ../lib . -name '*.hi' | xargs rm -f + find ../lib . -name '*.o' | xargs rm -f + find . -name '*.html' | xargs rm -f + rm -f test-mc test-mc.tix + rm -rf .hpc +