monte-carlo 0.1 → 0.2
raw patch · 17 files changed
+1363/−332 lines, 17 filesdep +arraydep +uvectordep ~gsl-random
Dependencies added: array, uvector
Dependency ranges changed: gsl-random
Files
- Control/Monad/MC.hs +0/−14
- Control/Monad/MC/GSL.hs +0/−271
- examples/Pi.hs +23/−18
- examples/Poker.hs +95/−0
- examples/Sampling.hs +58/−0
- lib/Control/Monad/MC.hs +14/−0
- lib/Control/Monad/MC/Base.hs +88/−0
- lib/Control/Monad/MC/Class.hs +34/−0
- lib/Control/Monad/MC/GSL.hs +34/−0
- lib/Control/Monad/MC/GSLBase.hs +285/−0
- lib/Control/Monad/MC/Repeat.hs +54/−0
- lib/Control/Monad/MC/Sample.hs +189/−0
- lib/Control/Monad/MC/Summary.hs +96/−0
- lib/Control/Monad/MC/Walker.hs +159/−0
- monte-carlo.cabal +45/−29
- tests/Main.hs +172/−0
- tests/Makefile +17/−0
− Control/Monad/MC.hs
@@ -1,14 +0,0 @@--------------------------------------------------------------------------------- |--- Module : Control.Monad.MC--- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>--- License : BSD3--- Maintainer : Patrick Perry <patperry@stanford.edu>--- Stability : experimental-----module Control.Monad.MC (- module Control.Monad.MC.GSL- ) where--import Control.Monad.MC.GSL
− Control/Monad/MC/GSL.hs
@@ -1,271 +0,0 @@-{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, UndecidableInstances #-}--------------------------------------------------------------------------------- |--- Module : Control.Monad.MC.GSL--- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>--- License : BSD3--- Maintainer : Patrick Perry <patperry@stanford.edu>--- 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)--}
examples/Pi.hs view
@@ -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
+ examples/Poker.hs view
@@ -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"
+ examples/Sampling.hs view
@@ -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)
+ lib/Control/Monad/MC.hs view
@@ -0,0 +1,14 @@+-----------------------------------------------------------------------------+-- |+-- Module : Control.Monad.MC+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module Control.Monad.MC (+ module Control.Monad.MC.GSL+ ) where++import Control.Monad.MC.GSL
+ lib/Control/Monad/MC/Base.hs view
@@ -0,0 +1,88 @@+{-# LANGUAGE TypeFamilies, MultiParamTypeClasses, FlexibleContexts #-}+-----------------------------------------------------------------------------+-- |+-- Module : Control.Monad.MC.Base+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- 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 #-}
+ lib/Control/Monad/MC/Class.hs view
@@ -0,0 +1,34 @@+-----------------------------------------------------------------------------+-- |+-- Module : Control.Monad.MC.Class+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- 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
+ lib/Control/Monad/MC/GSL.hs view
@@ -0,0 +1,34 @@+-----------------------------------------------------------------------------+-- |+-- Module : Control.Monad.MC.GSL+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- 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 )
+ lib/Control/Monad/MC/GSLBase.hs view
@@ -0,0 +1,285 @@+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, UndecidableInstances #-}+-----------------------------------------------------------------------------+-- |+-- Module : Control.Monad.MC.GSLBase+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- 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 #-}
+ lib/Control/Monad/MC/Repeat.hs view
@@ -0,0 +1,54 @@+-----------------------------------------------------------------------------+-- |+-- Module : Control.Monad.MC.Repeat+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- 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 #-}
+ lib/Control/Monad/MC/Sample.hs view
@@ -0,0 +1,189 @@+{-# LANGUAGE ScopedTypeVariables #-}+-----------------------------------------------------------------------------+-- |+-- Module : Control.Monad.MC.Sample+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- 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 #-}
+ lib/Control/Monad/MC/Summary.hs view
@@ -0,0 +1,96 @@+-----------------------------------------------------------------------------+-- |+-- Module : Control.Monad.MC.Summary+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- 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
+ lib/Control/Monad/MC/Walker.hs view
@@ -0,0 +1,159 @@+{-# LANGUAGE TypeOperators #-}+-----------------------------------------------------------------------------+-- |+-- Module : Control.Monad.MC.Walker+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- 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)+
monte-carlo.cabal view
@@ -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 <patperry@stanford.edu>-author: Patrick Perry-maintainer: Patrick Perry <patperry@stanford.edu>-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 <patperry@stanford.edu>+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
+ tests/Main.hs view
@@ -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)+ ]
+ tests/Makefile view
@@ -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+