monte-carlo 0.5 → 0.6

raw patch · 18 files changed

+714/−1006 lines, 18 filesdep +primitivedep +transformersdep −deepseqdep −mtldep ~QuickCheckdep ~gsl-randomPVP ok

version bump matches the API change (PVP)

Dependencies added: primitive, transformers

Dependencies removed: deepseq, mtl

Dependency ranges changed: QuickCheck, gsl-random

API changes (from Hackage documentation)

Files

• NEWS +32/−0
• examples/Binomial.hs +14/−12
• examples/Pi.lhs +22/−68
• examples/Poker.hs +15/−13
• examples/Queue.hs +82/−47
• lib/Control/Monad/MC.hs +1/−3
• lib/Control/Monad/MC/Base.hs +0/−181
• lib/Control/Monad/MC/Class.hs +0/−27
• lib/Control/Monad/MC/GSL.hs +6/−29
• lib/Control/Monad/MC/GSLBase.hs +155/−242
• lib/Control/Monad/MC/Repeat.hs +40/−17
• lib/Control/Monad/MC/Sample.hs +84/−199
• lib/Data/Summary.hs +0/−15
• lib/Data/Summary/Bool.hs +96/−50
• lib/Data/Summary/Double.hs +119/−57
• lib/Data/Summary/Utils.hs +9/−7
• monte-carlo.cabal +31/−32
• tests/Main.hs +8/−7

@@ -1,3 +1,35 @@+Changes in 0.6:++* Major overhaul. A lot of client code will break.++* Replaced the old MC type with a monad transformer,+  requiring the base monad to be an instance of PrimMonad.++* Removed the MCT type.++* Removed the MonadMC class.++* Removed unsafe operations relying on lazy IO (unsafeInterleaveIO),+  or replaced them by safe, strict versions.++* Added new fold operation foldMC.++* Refactored and simplified sampling and shuffling functions.  The+  "sampleSubset" functions are strict now.++* Hide Data.Summary.Utils, and remove Data.Summary.++* Change Data.Summary.Bool and Data.Summary.Double interfaces to mimic+  the Data.Set functions.  Like Data.Set, these modules should now+  be used with qualified imports.++* Removed NFData instances from Summary types; add Eq instances.++* Add Eq, Show, Data, Typeable instances.++* Updated examples.++ Changes in 0.5:  * Clark Gaebel added Monoid and NFData instances to Summary types

@@ -2,15 +2,15 @@     where  import Control.Monad+import Control.Monad.Primitive( PrimMonad ) import Data.List( foldl' ) import Text.Printf( printf )  import Control.Monad.MC-import Data.Summary-import Data.Summary.Utils( inInterval )+import qualified Data.Summary.Double as S  -- | Sample from a binomial distribution with the given parameters.-binomial :: (MonadMC m) => Int -> Double -> m Int+binomial :: (PrimMonad m) => Int -> Double -> MC m Int binomial n p = let     q     = 1 - p     probs = map (\i -> (fromIntegral $ n `choose` i) * p^^i * q^^(n-i)) [0..n]@@ -18,29 +18,31 @@  -- | 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 :: (PrimMonad m) => Int -> Double -> Int -> MC m (Double,Double) binomialMean n p reps =-    liftM (sampleCI 0.95 . summary . map fromIntegral) $-        replicateMC reps (binomial n p)+    liftM (S.meanCI 0.95) $+        foldMC (\s x -> return $! S.insertWith fromIntegral x s) S.empty+               reps (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 :: (PrimMonad m) => Int -> Double -> Int -> Int -> MC m Int coverage n p size reps =-    liftM (length . filter (mu `inInterval`)) $-        replicateMC reps $-            binomialMean n p size+    foldMC (\tot ci -> return $! update tot (mu `inInterval` ci)) 0+           reps (binomialMean n p size)   where     mu = fromIntegral n * p+    x `inInterval` (a,b) = x > a && x < b+    update tot b = tot + (if b then 1 else 0)  main =     let seed = 0-        reps = 100+        reps = 10000         n    = 10         p    = 0.2         size = 500-        c    =  coverage n p size reps `evalMC` mt19937 seed +        c    = (coverage n p size reps) `evalMC` (mt19937 seed)     in         printf "\nOf %d 95%%-intervals, %d contain the true value.\n" reps c 

@@ -4,11 +4,11 @@ of them fall in the unit circle.  \begin{code}-import Control.Monad( liftM, liftM2 )-import Control.Monad.MC( MC, uniform, replicateMC, evalMC, mt19937 )-import Data.Summary.Bool( summary, sampleMean, sampleSE )-import Data.Summary.Utils( interval )-import Text.Printf+import Control.Monad( liftM2 )+import Control.Monad.MC( STMC, evalMC, foldMC, mt19937, uniform )+import Data.Monoid( (<>), mempty )+import qualified Data.Summary.Bool as S+import Text.Printf( printf ) \end{code}  First, we need a function to test whether or not a point is in the@@ -19,75 +19,26 @@ inUnitCircle (x,y) = x*x + y*y <= 1 \end{code} -The first line is the type signature, which tells us that "inUnitCircle"-is a function which takes a pair of `Double`s and returns a `Bool`.  In-English, "::" means "has type".  The second line is the one that defines-the function.--\begin{code}-estimatePi :: [(Double,Double)] -> (Double,Double)-estimatePi xs =-    let s       = summary $ map inUnitCircle xs-        (mu,se) = (sampleMean s, sampleSE s) in-    (4*mu,4*se)-\end{code}- Next, we need to generate a random point in the unit box.  In the Control.Monad.MC module, there is a function for generating uniform-values in an interval, called "uniform".  This function has a +values in an interval, called "uniform".  This function has a funny-looking type, but you can think of it as: -    uniform :: Double -> Double -> MC Double+    uniform :: Double -> Double -> STMC s Double -This type means that the function takes the two endpoints of the +This type means that the function takes the two endpoints of the interval as arguments, and returns a Monte-Carlo action which produces a Double. -You can think of the type "MC Double" as a random number generator.-For general types, "MC a" is a generator for values of type "a".  In-fact, "MC" is one of a general class of objects called a Monad.  We-use Monads in Haskell for making sure events happen in the right order.-That is, if we have three parts of a simulation, say A, B, and C, and-we want them to happen in the order- -    A ====> B ====> C--then we would like to make sure that A is done consuming random numbers-before B consumes anything.  Likewise, we want B to finish before C-starts.  Monads are the magic that enable us to ensure this.--There are a number of functions in the standard library for working with-monads.  The first we will use is--    liftM2 :: (Monad m) => (a1 -> a2 -> r) -> m a1 -> m a2 -> m r--This function works on *any* Monad.  When we use it on the MC monad, it-will have type--    liftM2 :: (a1 -> a2 -> r) -> MC a1 -> MC a2 -> MC r-    -What liftM2 does is it takes a function of two arguments and two Monte-Carlo actions.  It returns a new Monte Carlo action that does the following:--  1. generate a random value of type a1 using the first action-  2. generate a random value of type a2 using the second action-  3. apply a function to the two values and return the result-  -We do not need to write liftM2 ourselves, since it is provided in the -"Control.Monad" module.  But, if we did have to define it, the code would-look like:+You can think of the type "STMC s Double" as a random number generator.+For general types, "STMC s a" is a generator for values of type "a". -liftM2 f ma1 ma2 = do-    a1 <- ma1-    a2 <- ma2-    return (f a1 a2)-    -This code for this uses the "do" notation of Haskell, which allows us-to specify a series of actions in sequential order.+The following code generates an x value in the range (-1,1), then a y+value in the range (-1,1), then returns the pair (x,y):  \begin{code}-unitBox :: MC (Double,Double)-unitBox = liftM2 (,) (uniform (-1) 1) +unitBox :: STMC s (Double,Double)+unitBox = liftM2 (,) (uniform (-1) 1)                      (uniform (-1) 1) \end{code} @@ -95,17 +46,20 @@ -- the sample mean and standard error.  \begin{code}-simulation :: Int -> MC (Double,Double)-simulation n = -    estimatePi `fmap` replicateMC n unitBox+estimatePi :: Int -> STMC s (Double,Double)+estimatePi n = let+    circ = fmap inUnitCircle unitBox+    mc   = foldMC (\s x -> return $! S.insert x s) S.empty n circ+    in fmap (\s -> let (mu,se) = (S.mean s, S.meanSE s)+                   in (4*mu,4*se)) mc \end{code}  \begin{code} main =     let seed    = 0         n       = 1000000 -        (mu,se) = simulation n `evalMC` mt19937 seed-        (l,u)   = interval 0.95 mu se+        (mu,se) = estimatePi n `evalMC` mt19937 seed+        (l,u)   = (mu - 2.576 * se, mu + 2.576 * se)     in do         printf "\nEstimate: %g" mu         printf "\n99%% Confidence Interval: (%g, %g)" l u

@@ -1,16 +1,18 @@ module Main where-    + import Control.Monad+import Control.Monad.Primitive( PrimMonad )  import Control.Monad.MC+import Control.Monad.ST import Data.List import Data.Map( Map ) import qualified Data.Map as Map 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 +data Card = Card { number :: Int                  , suit   :: Suit                  }           deriving (Eq, Show)@@ -24,19 +26,19 @@  -- | Get a list of cards that make up a 52-card deck. deck :: [Card]-deck = [ Card i s +deck = [ Card i s        | i <- [ ace..king ]        , s <- [ Club, Diamond, Heart, Spade ] ]  -- | A type for the various poker hands. data Hand = HighCard  | Pair | TwoPair | ThreeOfAKind | Straight | Flush-          | FullHouse | FourOfAKind | StraightFlush +          | 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 +hand cs =+    case matches of         [1,1,1,1,1] -> case undefined of                            _ | isStraight && isFlush -> StraightFlush                            _ | isFlush               -> Flush@@ -50,7 +52,7 @@   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 ] @@ -58,9 +60,9 @@      matches = (sort . map length . group) (x:xs) -    + -- | Deal a five-card hand by choosing a random subset of the deck.-deal :: (MonadMC m) => m [Card]+deal :: (PrimMonad m) => MC m [Card] deal = sampleSubset deck 5  -- | A type for storing the frequencies of the various hands.@@ -78,15 +80,15 @@ main =     let seed   = 0         reps   = 100000-        counts = foldl' updateCounts emptyCounts $ -                     replicateMC reps deal `evalMC` mt19937 seed +        counts = foldl' updateCounts emptyCounts $+                     replicateMC reps deal (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 +                p     = fromIntegral c / n                 se    = sqrt (p * (1 - p) / n)                 delta = 2.575829 * se                 (l,u) = (p-delta, p+delta) in

@@ -1,34 +1,44 @@  import Control.Monad-import Control.Monad.MC+import Control.Monad.Primitive( PrimMonad ) import Data.List( foldl' )-import Data.Summary import Text.Printf( printf ) +import Control.Monad.MC+import Data.Summary.Double( Summary )+import qualified Data.Summary.Double as S++ -- | There a three items on the menu. data Item = Cheeseburger | Fries | Milkshake + -- | A customer orders some number of items data Customer = Customer { orderOf :: [Item] } + -- | The order size is a Poisson random variable with mean 2.-orderSize :: MC Int+orderSize :: (PrimMonad m) => MC m Int orderSize = liftM (1+) $ poisson 2 + -- | The items are sampled with the given weights.-item :: MC Item+item :: (PrimMonad m) => MC m Item item = sampleWithWeights [ (4, Cheeseburger), (2, Fries), (1, Milkshake) ] + -- | Generate a random order.-order :: MC [Item]+order :: (PrimMonad m) => MC m [Item] order = do     n <- orderSize     replicateM n item + -- | Generate a random customer.-customer :: MC Customer+customer :: (PrimMonad m) => MC m Customer customer = liftM Customer order-    ++ -- | A customer event.  The interarrival time is the time that elapeses -- between when the previous customer arrives and when the current customer  -- arrives.@@ -36,32 +46,37 @@                                    , interarrivalTime :: !Double                                    } + -- | Generate a random customer event.  The interarrival time distribution -- is exponential with mean 1.-customerEvent :: MC CustomerEvent+customerEvent :: (PrimMonad m) => MC m CustomerEvent customerEvent = do     c     <- customer     delta <- exponential 10     return $ CustomerEvent c delta + -- | The time it takes to make an item.-cook :: Item -> MC Double+cook :: (PrimMonad m) => Item -> MC m Double cook Cheeseburger = exponential 3 cook Fries        = exponential 1 cook Milkshake    = exponential 2 + -- | The time it takes to cook all of the items in the list is equal -- to the maximum time.-cookAll :: [Item] -> MC Double+cookAll :: (PrimMonad m) => [Item] -> MC m Double cookAll items = do     ts <- mapM cook items     return $ foldl' max 0 ts + -- | A customer in line, along with how long they have been waiting. data Waiting = Waiting { waiting        :: !Customer                        , hasBeenWaiting :: !Double                        } + -- | A customer, along with how long it takes to prepare the customer's order -- and how long the customer has to wait. data Service = Service { serving     :: !Customer@@ -69,37 +84,49 @@                        , serviceTime :: !Double                        } + -- | Given a customer who has been wating in line, provide them with service. -- If the customer has been waiting for longer than 5 minutes, work twice as -- fast to cook the food.-serveWaiting :: Waiting -> MC Service+serveWaiting :: (PrimMonad m) => Waiting -> MC m Service serveWaiting (Waiting c w) = do     t <- cookAll $ orderOf c     let t' = if w > 5 then 0.5*t else t     return $ Service c w t' + -- | A resturant has one server, who may be busy. There is a list of -- customers wating in line. data Restaurant = Restaurant { inProgress  :: Maybe InProgress                              , waitingLine :: [Waiting]                              } + -- | An in-progress service event. data InProgress = InProgress { service      :: !Service                              , timeToFinish :: !Double                              } ++-- | An empty restaurant.+emptyRestaurant :: Restaurant+emptyRestaurant = Restaurant Nothing []++ -- | Update the amount of time the customers have been waiting by adding -- the given amount. addToWait :: Double -> [Waiting] -> [Waiting] addToWait delta = map (\(Waiting w t) -> Waiting w (t+delta)) + -- | Serve customers in the restaurant for the given amount of time.        -serveForTime :: Double -> Restaurant -> MC ([Service], Restaurant)-serveForTime = +serveForTime :: (PrimMonad m) => Double+                              -> Restaurant+                              -> MC m ([Service], Restaurant)+serveForTime =     let serveForTimeHelp ss t r = case r of             -- When no one is being served and no one is in line, do nothing.-            Restaurant Nothing  [] -> +            Restaurant Nothing  [] ->                 return $ (ss, r)              -- When no one is being served, take the first person in line@@ -124,66 +151,74 @@                                        xs'    = addToWait t xs                                        r'     = Restaurant y' xs' in                                    return (ss,r')-    in serveForTimeHelp []                                   +    in serveForTimeHelp [] + -- | Given a new customer arrival event, produce a list of all of the new -- service events that happen before the customer gets there, and return -- the updated restaurant state at the time immediately after the customer -- arrives.-processEvent :: CustomerEvent    -             -> Restaurant       -             -> MC ([Service], Restaurant)+processEvent :: (PrimMonad m) => CustomerEvent+                              -> Restaurant+                              -> MC m ([Service], Restaurant) processEvent (CustomerEvent c t) r = do     (ss,(Restaurant y xs)) <- serveForTime t r     return $ (ss, (Restaurant y $ xs ++ [Waiting c 0])) + -- | Finish serving all of the customers in line.-finishServing :: Restaurant -> MC [Service]+finishServing :: (PrimMonad m) => Restaurant -> MC m [Service] finishServing r = do     (ss,_) <- serveForTime infinity r     return ss   where     infinity = 1/0-    --- | A restaurant takes a list of customer events and generates a random--- list of service events.  The reason for the call to "unsafeInterleaveMC"--- is that we want to make sure that we return a lazy list.   Without it,--- the function will return only after it has consumed all of the random--- numbers it needs.  This is problemeatic if the input list is large or--- or infinite.-restaurant :: [CustomerEvent] -> MC [Service]-restaurant = -    let restaurantHelp r []     = finishServing r-        restaurantHelp r (c:cs) = unsafeInterleaveMC $ do-            (ss,r') <- processEvent c r-            ss'     <- restaurantHelp r' cs-            return $ ss ++ ss'-    in restaurantHelp (Restaurant Nothing []) ++-- | Run a restaurnt.  Whenever a new set of service events is generated,+-- update the accumulator.+foldRestaurant :: (PrimMonad m) => (a -> [Service] -> MC m a)+                                -> a+                                -> [CustomerEvent]+                                -> Restaurant+                                -> MC m a+foldRestaurant f a []     r = finishServing r >>= f a+foldRestaurant f a (c:cs) r = do+    (ss,r') <- processEvent c r+    a' <- f a ss+    foldRestaurant f a' cs r'+++-- | Compute a summary of the total waiting times for each customer.+summarizeService :: (PrimMonad m) => [CustomerEvent] -> Restaurant -> MC m Summary+summarizeService cs r =+    foldRestaurant (\s ss -> return $!+                                 foldl' (flip $ S.insertWith totalTime) s ss)+                   S.empty cs r+  where+    totalTime (Service _ w s) = w+s++ -- | An infinite stream of customerEvents.  This stream uses its own private  -- random number generator (mt19937 is the Mersenne-Twister algorithm). customerEvents :: Seed -> [CustomerEvent]-customerEvents seed = repeatMC customerEvent `evalMC` mt19937 seed+customerEvents seed = customerEvent `repeatMC` (mt19937 seed) + -- | Given a seed for the customers and a seed for the restaurant, run the -- simulation.-simulation :: Seed -> Seed -> [Service]-simulation customerSeed restaurantSeed =-    restaurant (customerEvents customerSeed) `evalMC` mt19937 restaurantSeed+simulation :: Seed -> Seed -> Int -> Summary+simulation customerSeed restaurantSeed n = let+    cs = take n $ customerEvents customerSeed+    r  = emptyRestaurant+    in evalMC (summarizeService cs r) (mt19937 restaurantSeed) --- | Compute a summary of the total waitings time for each customer.-summarize :: [Service] -> Summary-summarize = summary . map totalTime-  where -    totalTime (Service _ w s) = w+s-     -- | Run the program-main = +main =     let customerSeed    = 0         restaurantSeed  = 100         numTransactions = 100000-        results         = summarize $ take numTransactions $ -                              simulation  customerSeed restaurantSeed+        results         = simulation customerSeed restaurantSeed numTransactions     in do         putStrLn ""         putStrLn "Total Service Time:"

@@ -6,9 +6,7 @@ -- Maintainer : Patrick Perry <patperry@gmail.com> -- Stability  : experimental ----- A monad and monad transformer for monte carlo computations.  Currently,--- the default is the GNU Scientific Library-based implementation, but this--- may change in the future.+-- A monad and monad transformer for Monte Carlo computations. --  module Control.Monad.MC (

@@ -1,181 +0,0 @@-{-# LANGUAGE TypeFamilies, PolyKinds #-}--------------------------------------------------------------------------------- |--- Module     : Control.Monad.MC.Base--- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>--- License    : BSD3--- Maintainer : Patrick Perry <patperry@gmail.com>--- Stability  : experimental-----module Control.Monad.MC.Base-    where--import Control.Monad-import qualified Control.Monad.MC.GSLBase as GSL--import qualified Data.Vector.Storable as VS--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--    -- | @exponential mu@ generates an Exponential variate with mean @mu@.-    exponential :: Double -> m Double--    -- | @levy c alpha@ gets a Levy alpha-stable variate with scale @c@ and-    -- exponent @alpha@.  The algorithm only works for @0 < alpha <= 2@.-    levy :: Double -> Double -> m Double--    -- | @levySkew c alpha beta @ gets a skew Levy alpha-stable variate-    -- with scale @c@, exponent @alpha@, and skewness @beta@.  The skew-    -- parameter must lie in the range @[-1,1]@.  The algorithm only works-    -- for @0 < alpha <= 2@.-    levySkew :: Double -> Double -> Double -> m Double--    -- | @poisson mu@ generates a Poisson random variable with mean @mu@.-    poisson :: Double -> m Int--    -- | @cauchy a@ generates a Cauchy random variable with scale-    -- parameter @a@.-    cauchy :: Double -> m Double--    -- | @beta a b@ generates a beta random variable with-    -- parameters @a@ and @b@.-    beta :: Double -> Double -> m Double--    -- | @logistic a@ generates a logistic random variable with-    -- parameter @a@.-    logistic :: Double -> m Double--    -- | @pareto a b@ generates a Pareto random variable with-    -- exponent @a@ and scale @b@.-    pareto :: Double -> Double -> m Double--    -- | @weibull a b@ generates a Weibull random variable with-    -- scale @a@ and exponent @b@.-    weibull :: Double -> Double -> m Double--    -- | @gamma a b@ generates a gamma random variable with-    -- parameters @a@ and @b@.-    gamma :: Double -> Double -> m Double--    -- | @multinomial n ps@ generates a multinomial random-    -- variable with parameters @ps@ formed by @n@ trials.-    multinomial :: Int -> VS.Vector Double -> m (VS.Vector Int)--    -- | @dirichlet alphas@ generates a Dirichlet random variable-    -- with parameters @alphas@.-    dirichlet :: VS.Vector Double -> m (VS.Vector Double)--    -- | Get the baton from the Monte Carlo monad without performing any-    -- computations.  Useful but dangerous.-    unsafeInterleaveMC :: m a -> m a----- | Generate 'True' events with the given probability-bernoulli :: (MonadMC m) => Double -> m Bool-bernoulli p = liftM (< p) $ uniform 0 1-{-# INLINE bernoulli #-}--------------------------------- 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 #-}-    exponential = GSL.exponential-    {-# INLINE exponential #-}-    levy = GSL.levy-    {-# INLINE levy #-}-    levySkew = GSL.levySkew-    {-# INLINE levySkew #-}-    poisson = GSL.poisson-    {-# INLINE poisson #-}-    cauchy = GSL.cauchy-    {-# INLINE cauchy #-}-    beta = GSL.beta-    {-# INLINE beta #-}-    logistic = GSL.logistic-    {-# INLINE logistic #-}-    pareto = GSL.pareto-    {-# INLINE pareto #-}-    weibull = GSL.weibull-    {-# INLINE weibull #-}-    gamma = GSL.gamma-    {-# INLINE gamma #-}-    multinomial = GSL.multinomial-    {-# INLINE multinomial #-}-    dirichlet = GSL.dirichlet-    {-# INLINE dirichlet #-}-    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 #-}-    exponential mu = GSL.liftMCT $ GSL.exponential mu-    {-# INLINE exponential #-}-    levy c alpha = GSL.liftMCT $ GSL.levy c alpha-    {-# INLINE levy #-}-    levySkew c alpha beta_ = GSL.liftMCT $ GSL.levySkew c alpha beta_-    {-# INLINE levySkew #-}-    poisson mu = GSL.liftMCT $ GSL.poisson mu-    {-# INLINE poisson #-}-    cauchy a = GSL.liftMCT $ GSL.cauchy a-    {-# INLINE cauchy #-}-    beta a b = GSL.liftMCT $ GSL.beta a b-    {-# INLINE beta #-}-    logistic a = GSL.liftMCT $ GSL.logistic a-    {-# INLINE logistic #-}-    pareto a b = GSL.liftMCT $ GSL.pareto a b-    {-# INLINE pareto #-}-    weibull a b = GSL.liftMCT $ GSL.weibull a b-    {-# INLINE weibull #-}-    gamma a b = GSL.liftMCT $ GSL.gamma a b-    {-# INLINE gamma #-}-    multinomial n ps = GSL.liftMCT $ GSL.multinomial n ps-    {-# INLINE multinomial #-}-    dirichlet alphas = GSL.liftMCT $ GSL.dirichlet alphas-    {-# INLINE dirichlet #-}-    unsafeInterleaveMC = GSL.unsafeInterleaveMCT-    {-# INLINE unsafeInterleaveMC #-}

@@ -1,27 +0,0 @@--------------------------------------------------------------------------------- |--- Module     : Control.Monad.MC.Class--- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>--- License    : BSD3--- Maintainer : Patrick Perry <patperry@gmail.com>--- Stability  : experimental------ The abstract MonadMC interface and utility functions for Monte Carlo--- computations.-----module Control.Monad.MC.Class (-    -- * The Monte Carlo monad type class-    HasRNG(..),-    MonadMC(..),--    -- * Random distributions-    bernoulli,--    module Control.Monad.MC.Sample,-    module Control.Monad.MC.Repeat,-    ) where--import Control.Monad.MC.Base-import Control.Monad.MC.Sample-import Control.Monad.MC.Repeat

@@ -10,34 +10,11 @@ -- of the functions in the GNU Scientific Library.  module Control.Monad.MC.GSL (-    -- * The Monte Carlo monad-    MC,-    runMC,-    evalMC,-    execMC,--    -- * The Monte Carlo monad transformer-    MCT,-    runMCT,-    evalMCT,-    execMCT,-    liftMCT,--    -- * Pure random number generator creation-    RNG,-    Seed,-    mt19937,-    mt19937WithState,-    rngName,-    rngSize,-    rngState,--    -- * Overloaded Monte Carlo monad interface-    module Control.Monad.MC.Class,-+    module Control.Monad.MC.GSLBase,+    module Control.Monad.MC.Sample,+    module Control.Monad.MC.Repeat,     ) where -import Control.Monad.MC.GSLBase ( MC, runMC, evalMC, execMC,-    MCT, runMCT, evalMCT, execMCT, liftMCT, RNG, Seed, mt19937, mt19937WithState,-    rngName, rngSize, rngState )-import Control.Monad.MC.Class hiding ( RNG )+import Control.Monad.MC.GSLBase+import Control.Monad.MC.Sample+import Control.Monad.MC.Repeat

@@ -1,4 +1,4 @@-{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses, UndecidableInstances #-}+{-# LANGUAGE DeriveDataTypeable, RankNTypes, TypeFamilies #-} ----------------------------------------------------------------------------- -- | -- Module     : Control.Monad.MC.GSLBase@@ -9,99 +9,99 @@ --  module Control.Monad.MC.GSLBase (-    -- * The Monte Carlo monad+    -- * Monte Carlo monad transformer     MC(..),-    runMC,+    STMC,+    IOMC,     evalMC,-    execMC,-    unsafeInterleaveMC, -    -- * The Monte Carlo monad transformer-    MCT(..),-    runMCT,-    evalMCT,-    execMCT,-    unsafeInterleaveMCT,-    liftMCT,--    -- * Pure random number generator creation+    -- * Random number generator+    -- ** Types     RNG,+    IORNG,+    STRNG,     Seed,+    -- ** Creation     mt19937,     mt19937WithState,-    rngName,-    rngSize,-    rngState,--    -- * Getting and setting the random number generator-    getRNG,-    setRNG,+    -- ** State+    getRNGName,+    getRNGSize,+    getRNGState,+    setRNGState,      -- * Random number distributions+    -- ** Uniform     uniform,     uniformInt,+    -- ** Continuous     normal,     exponential,+    gamma,+    cauchy,     levy,     levySkew,-    poisson,-    cauchy,-    beta,-    logistic,     pareto,     weibull,-    gamma,-    multinomial,+    logistic,+    beta,+    -- ** Discrete+    bernoulli,+    poisson,+    -- ** Multivariate     dirichlet,+    multinomial,     ) where -import Control.Applicative      ( Applicative(..), (<$>) )-import Control.Monad            ( ap, 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, unsafeInterleaveIO )+import Control.Applicative       ( Applicative(..) )+import Control.Monad             ( liftM )+import Control.Monad.Fix         ( MonadFix(..) )+import Control.Monad.IO.Class    ( MonadIO(..) )+import Control.Monad.ST          ( ST, runST )+import Control.Monad.Primitive   ( PrimMonad(..), unsafePrimToPrim )+import Control.Monad.Trans.Class ( MonadTrans(..) )+import Data.Typeable             ( Typeable )+import Data.Word                 ( Word8, Word64 )  import qualified Data.Vector.Storable as VS  import qualified GSL.Random.Gen as GSL import GSL.Random.Dist --- | A Monte Carlo monad with an internal random number generator.-newtype MC a = MC (GSL.RNG -> IO a) --- | 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) (RNG r) = unsafePerformIO $ do-    r' <- GSL.cloneRNG r-    a  <- g r'-    return (a,RNG r')-{-# NOINLINE runMC #-}+-- | A Monte Carlo monad transformer.  This type provides access+-- to a random number generator while allowing operations in a+-- base monad, @m@.+newtype MC m a = MC { runMC :: RNG (PrimState m) -> m a } --- | 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+-- | Type alias for when the base monad is 'ST'.+type STMC s a = MC (ST s) a -unsafeInterleaveMC :: MC a -> MC a-unsafeInterleaveMC (MC m) = MC $ \r ->-    unsafeInterleaveIO (m r)+-- | Type alias for when the base monad is 'IO'.+type IOMC a = MC IO a -instance Functor MC where-    fmap f (MC m) = MC $ \r ->-        fmap f (m r) -instance Monad MC where+-- | Evaluate the result of a Monte Carlo computation using the given+-- random number generator.+evalMC :: (forall s. STMC s a) -> (forall s. ST s (STRNG s)) -> a+evalMC ma mr = runST $ do+    r <- mr+    runMC ma r+++instance (Functor m) => Functor (MC m) where+    fmap f (MC m) = MC $ \r -> fmap f (m r)+    {-# INLINE fmap #-}++instance (Applicative m) => Applicative (MC m) where+    pure a = MC $ \_ -> pure a+    {-# INLINE pure #-}++    (MC gf) <*> (MC ga) = MC $ \r -> (gf r) <*> (ga r)+    {-# INLINE (<*>) #-}++instance (Monad m) => Monad (MC m) where     return a = MC $ \_ -> return a     {-# INLINE return #-} @@ -111,250 +111,163 @@             in m' r     {-# INLINE (>>=) #-} -    fail s = MC $ \_ -> fail s-    {-# INLINE fail #-}--instance Applicative MC where-    pure  = return-    (<*>) = ap---- | A parameterizable Monte Carlo monad for encapsulating an inner--- monad.-newtype MCT m a = MCT (GSL.RNG -> IO (m a))---- | Similar to 'runMC'.-runMCT :: (Monad m) => MCT m a -> RNG -> m (a,RNG)-runMCT (MCT g) (RNG r) = unsafePerformIO $ do-    r' <- GSL.cloneRNG r-    ma <- g r'-    return (ma >>= \a -> return (a, RNG 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 g) = MCT $ \r -> do-    a <- g r-    return (return a)-{-# INLINE liftMCT #-}--unsafeInterleaveMCT :: (Monad m) => MCT m a -> MCT m a-unsafeInterleaveMCT (MCT g) = MCT $ \r ->-    unsafeInterleaveIO (g r)-{-# INLINE unsafeInterleaveMCT #-}--instance (Monad m) => Functor (MCT m) where-    fmap f (MCT g) = MCT $ \r -> do-        ma <- g r-        return (ma >>= return . f)-    {-# INLINE fmap #-}--instance (Monad m) => Monad (MCT m) where-    return a = MCT $ \_ -> return (return a)-    {-# INLINE return #-}--    (MCT g) >>= k =-        MCT $ \r -> do-            ma <- g r-            return $ ma >>= \a ->-                let (MCT m') = k a-                in unsafePerformIO $ m' r-    {-# NOINLINE (>>=) #-}--    fail str = MCT $ \_ -> fail str+    fail msg = MC $ \_ -> fail msg     {-# INLINE fail #-} -instance (Monad m) => Applicative (MCT m) where-    pure  = return-    (<*>) = ap--instance (MonadPlus m) => MonadPlus (MCT m) where-    mzero = MCT $ \_ -> mzero-    {-# INLINE mzero #-}+instance (MonadFix m) => MonadFix (MC m) where+    mfix f = MC $ \r -> mfix $ flip (runMC . f) r+    {-# INLINE mfix #-} -    (MCT m) `mplus` (MCT n) =-        MCT $ \r -> do-            r' <- GSL.cloneRNG r-            mr <- m r-            nr <- n r'-            return (mr `mplus` nr)+instance (MonadIO m) => MonadIO (MC m) where+    liftIO io = MC $ \_ -> liftIO io+    {-# INLINE liftIO #-} -instance MonadTrans MCT where-    lift m = MCT $ \_ -> return m+instance MonadTrans MC where+    lift m = MC $ \_ -> m     {-# INLINE lift #-} -instance (MonadCont m) => MonadCont (MCT m) where-    callCC f = MCT $ \r ->-        return $ callCC $ \k ->-            let (MCT m) = f (\a -> MCT $ \_ -> return (k a))-            in unsafePerformIO (m r)-    {-# NOINLINE callCC #-} -instance (MonadError e m) => MonadError e (MCT m) where-    throwError             = lift . throwError-    {-# INLINE throwError #-}--    (MCT g) `catchError` h = MCT $ \r -> do-        ma <- g r-        return $ ma `catchError` \e ->-            let (MCT m') = h e-            in unsafePerformIO (m' r)-    {-# NOINLINE catchError #-}--instance (MonadIO m) => MonadIO (MCT m) where-    liftIO = lift . liftIO-    {-# INLINE liftIO #-}--instance (MonadReader r m) => MonadReader r (MCT m) where-    ask              = lift ask-    {-# INLINE ask #-}--    local f (MCT g) = MCT $ \r -> do-        ma <- g r-        return $ local f ma-    {-# INLINE local #-}--instance (MonadState s m) => MonadState s (MCT m) where-    get = lift get-    {-# INLINE get #-}--    put = lift . put-    {-# INLINE put #-}--instance (MonadWriter w m) => MonadWriter w (MCT m) where-    tell           = lift . tell-    {-# INLINE tell #-}--    listen (MCT g) = MCT $ \r -> do-        ma <- g r-        return (listen ma)-    {-# INLINE listen #-}+---------------------------- Random Number Generators ----------------------- -    pass (MCT g) = MCT $ \r -> do-        maf <- g r-        return (pass maf)-    {-# INLINE pass #-}+-- | The random number generator type.+newtype RNG s = RNG GSL.RNG+    deriving(Eq, Show, Typeable) ----------------------------- Random Number Generators -----------------------+-- | A shorter name for RNG in the 'IO' monad.+type IORNG = RNG (PrimState IO) --- | The random number generator type associated with 'MC' and 'MCT'.-newtype RNG = RNG GSL.RNG+-- | A shorter name for RNG in the 'ST' monad.+type STRNG s = RNG (PrimState (ST s))  -- | The seed type for the random number generators. type Seed = Word64 + -- | Get the name of the random number generator algorithm.-rngName :: RNG -> String-rngName (RNG r) = unsafePerformIO $ GSL.getName r-{-# NOINLINE rngName #-}+getRNGName :: (PrimMonad m) => RNG (PrimState m) -> m String+getRNGName (RNG r) = unsafePrimToPrim $ GSL.getName r  -- | Get the size of the generator state, in bytes.-rngSize :: RNG -> Int-rngSize (RNG r) = fromIntegral $ unsafePerformIO $ GSL.getSize r-{-# NOINLINE rngSize #-}+getRNGSize :: (PrimMonad m) => RNG (PrimState m) -> m Int+getRNGSize (RNG r) = liftM fromIntegral $ unsafePrimToPrim $ GSL.getSize r  -- | Get the state of the generator.-rngState :: RNG -> [Word8]-rngState (RNG r) = unsafePerformIO $ GSL.getState r-{-# NOINLINE rngState #-}--getRNG :: MC RNG-getRNG = MC (\r -> liftM RNG $ GSL.cloneRNG r)-{-# INLINE getRNG #-}+getRNGState :: (PrimMonad m) => RNG (PrimState m) -> m [Word8]+getRNGState (RNG r) = unsafePrimToPrim $ GSL.getState r -setRNG :: RNG -> MC ()-setRNG (RNG r') = MC $ \r -> GSL.copyRNG r r'-{-# INLINE setRNG #-}+-- | Set the state of the generator.+setRNGState :: (PrimMonad m) => RNG (PrimState m) -> [Word8] -> m ()+setRNGState (RNG r) xs = unsafePrimToPrim $ GSL.setState r xs --- | Get a Mersenne Twister random number generator seeded with the given+-- | Create a Mersenne Twister random number generator seeded with the given -- value.-mt19937 :: Seed -> RNG-mt19937 s = unsafePerformIO $ do+mt19937 :: (PrimMonad m) => Seed -> m (RNG (PrimState m))+mt19937 s = unsafePrimToPrim $ do     r <- GSL.newRNG GSL.mt19937     GSL.setSeed r s     return (RNG r)-{-# NOINLINE mt19937 #-} --- | Get a Mersenne Twister seeded with the given state.-mt19937WithState :: [Word8] -> RNG-mt19937WithState xs = unsafePerformIO $ do+-- | Create a Mersenne Twister seeded with the given state.+mt19937WithState :: (PrimMonad m) => [Word8] -> m (RNG (PrimState m))+mt19937WithState xs = unsafePrimToPrim $ do     r <- GSL.newRNG GSL.mt19937     GSL.setState r xs     return (RNG r)-{-# NOINLINE mt19937WithState #-}  -------------------------- Random Number Distributions ---------------------- -uniform :: Double -> Double -> MC Double+-- | @uniform a b@ generates a value uniformly distributed in @[a,b)@.+uniform :: (PrimMonad m) => Double -> Double -> MC m Double uniform 0 1 = liftRan0 GSL.getUniform uniform a b = liftRan2 getFlat a b -uniformInt :: Int -> MC Int+-- | @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 :: (PrimMonad m) => Int -> MC m Int uniformInt = liftRan1 GSL.getUniformInt -normal :: Double -> Double -> MC Double-normal 0  1     =            liftRan0 getUGaussianRatioMethod-normal mu 1     = (mu +) <$> liftRan0 getUGaussianRatioMethod-normal 0  sigma =            liftRan1 getGaussianRatioMethod sigma-normal mu sigma = (mu +) <$> liftRan1 getGaussianRatioMethod sigma+-- | @normal mu sigma@ generates a Normal random variable with mean+-- @mu@ and standard deviation @sigma@.+normal :: (PrimMonad m) => Double -> Double -> MC m Double+normal 0  1     =                liftRan0 getUGaussianRatioMethod+normal mu 1     = liftM (mu +) $ liftRan0 getUGaussianRatioMethod+normal 0  sigma =                liftRan1 getGaussianRatioMethod sigma+normal mu sigma = liftM (mu +) $ liftRan1 getGaussianRatioMethod sigma -exponential :: Double -> MC Double+-- | @exponential mu@ generates an Exponential variate with mean @mu@.+exponential :: (PrimMonad m) => Double -> MC m Double exponential = liftRan1 getExponential -poisson :: Double -> MC Int+-- | @poisson mu@ generates a Poisson random variable with mean @mu@.+poisson :: (PrimMonad m) => Double -> MC m Int poisson = liftRan1 getPoisson -levy :: Double -> Double -> MC Double+-- | @levy c alpha@ gets a Levy alpha-stable variate with scale @c@ and+-- exponent @alpha@.  The algorithm only works for @0 < alpha <= 2@.+levy :: (PrimMonad m) => Double -> Double -> MC m Double levy = liftRan2 getLevy -levySkew :: Double -> Double -> Double -> MC Double+-- | @levySkew c alpha beta @ gets a skew Levy alpha-stable variate+-- with scale @c@, exponent @alpha@, and skewness @beta@.  The skew+-- parameter must lie in the range @[-1,1]@.  The algorithm only works+-- for @0 < alpha <= 2@.+levySkew :: (PrimMonad m) => Double -> Double -> Double -> MC m Double levySkew = liftRan3 getLevySkew -cauchy :: Double -> MC Double+-- | @cauchy a@ generates a Cauchy random variable with scale+-- parameter @a@.+cauchy :: (PrimMonad m) => Double -> MC m Double cauchy = liftRan1 getCauchy -beta :: Double -> Double -> MC Double+-- | @beta a b@ generates a beta random variable with+-- parameters @a@ and @b@.+beta :: (PrimMonad m) => Double -> Double -> MC m Double beta = liftRan2 getBeta -logistic :: Double -> MC Double+-- | @logistic a@ generates a logistic random variable with+-- parameter @a@.+logistic :: (PrimMonad m) => Double -> MC m Double logistic = liftRan1 getLogistic -pareto :: Double -> Double -> MC Double+-- | @pareto a b@ generates a Pareto random variable with+-- exponent @a@ and scale @b@.+pareto :: (PrimMonad m) => Double -> Double -> MC m Double pareto = liftRan2 getPareto -weibull :: Double -> Double -> MC Double+-- | @weibull a b@ generates a Weibull random variable with+-- scale @a@ and exponent @b@.+weibull :: (PrimMonad m) => Double -> Double -> MC m Double weibull = liftRan2 getWeibull -gamma :: Double -> Double -> MC Double+-- | @gamma a b@ generates a gamma random variable with+-- parameters @a@ and @b@.+gamma :: (PrimMonad m) => Double -> Double -> MC m Double gamma = liftRan2 getGamma -multinomial :: Int -> VS.Vector Double -> MC (VS.Vector Int)+-- | @multinomial n ps@ generates a multinomial random+-- variable with parameters @ps@ formed by @n@ trials.+multinomial :: (PrimMonad m) => Int -> VS.Vector Double -> MC m (VS.Vector Int) multinomial = liftRan2 getMultinomial -dirichlet :: VS.Vector Double -> MC (VS.Vector Double)+-- | @dirichlet alphas@ generates a Dirichlet random variable+-- with parameters @alphas@.+dirichlet :: (PrimMonad m) => VS.Vector Double -> MC m (VS.Vector Double) dirichlet = liftRan1 getDirichlet +-- | Generate 'True' events with the given probability.+bernoulli :: (PrimMonad m ) => Double -> MC m Bool+bernoulli p = liftM (< p) $ uniform 0 1+{-# INLINE bernoulli #-}  -liftRan0 :: (GSL.RNG -> IO a) -> MC a-liftRan0 = MC+liftRan0 :: (PrimMonad m) => (GSL.RNG -> IO a) -> MC m a+liftRan0 f = MC $ \(RNG r) -> unsafePrimToPrim $ f r -liftRan1 :: (GSL.RNG -> a -> IO b) -> a -> MC b-liftRan1 f a = MC $ \r -> f r a+liftRan1 :: (PrimMonad m) => (GSL.RNG -> a -> IO b) -> a -> MC m b+liftRan1 f a = MC $ \(RNG r) -> unsafePrimToPrim $ f r a -liftRan2 :: (GSL.RNG -> a -> b -> IO c) -> a -> b -> MC c-liftRan2 f a b = MC $ \r -> f r a b+liftRan2 :: (PrimMonad m) => (GSL.RNG -> a -> b -> IO c) -> a -> b -> MC m c+liftRan2 f a b = MC $ \(RNG r) -> unsafePrimToPrim $ f r a b -liftRan3 :: (GSL.RNG -> a -> b -> c -> IO d) -> a -> b -> c -> MC d-liftRan3 f a b c = MC $ \r -> f r a b c+liftRan3 :: (PrimMonad m) => (GSL.RNG -> a -> b -> c -> IO d) -> a -> b -> c -> MC m d+liftRan3 f a b c = MC $ \(RNG r) -> unsafePrimToPrim $ f r a b c

@@ -1,3 +1,4 @@+{-# LANGUAGE RankNTypes #-} ----------------------------------------------------------------------------- -- | -- Module     : Control.Monad.MC.Repeat@@ -9,28 +10,50 @@  module Control.Monad.MC.Repeat (     -- * Repeating computations+    foldMC,     repeatMC,     replicateMC,     ) where -import Control.Monad.MC.Base+import Control.Monad.Primitive( PrimMonad )+import Control.Monad.MC.GSLBase+import Control.Monad.ST( ST, runST )+import Control.Monad.ST.Unsafe( unsafeInterleaveST ) --- | Produce a lazy infinite list of values from the given Monte Carlo--- generator.-repeatMC :: (MonadMC m) => m a -> m [a]-repeatMC = interleaveSequence . repeat-{-# INLINE repeatMC #-} +-- | Generate a sequence of replicates and incrementally consume+-- them via a left fold.+--+-- This fold is /not/ strict. The replicate consumer is responsible for+-- forcing the evaluation of its result to avoid space leaks.+foldMC :: (PrimMonad m) => (a -> b -> MC m a) -- ^ Replicate consumer.+                        -> a                  -- ^ Initial state for replicate consumer.+                        -> Int                -- ^ Number of replicates.+                        -> MC m b             -- ^ Generator.+                        -> MC m a+foldMC f a n mb | n <= 0    = return a+                | otherwise = do+    b <- mb+    a' <- f a b+    foldMC f a' (n-1) mb+{-# INLINE foldMC #-}+++-- | Produce a lazy infinite list of replicates from the given random+-- number generator and Monte Carlo procedure.+repeatMC :: (forall s. STMC s a) -> (forall s. ST s (STRNG s)) -> [a]+repeatMC mc mrng = runST $ do+    rng <- mrng+    go $ runMC mc rng+  where+    go m = unsafeInterleaveST $ do+        a  <- m+        as <- go m+        return (a:as)++ -- | Produce a lazy list of the given length using the specified--- generator.-replicateMC :: (MonadMC m) => Int -> m a -> m [a]-replicateMC n = interleaveSequence . replicate n-{-# INLINE replicateMC #-}+-- random number genrator and Monte Carlo procedure.+replicateMC :: Int -> (forall s. STMC s a) -> (forall s. ST s (STRNG s)) -> [a]+replicateMC n mc mrng = take n $ repeatMC mc mrng -interleaveSequence :: (MonadMC m) => [m a] -> m [a]-interleaveSequence []     = return []-interleaveSequence (m:ms) = unsafeInterleaveMC $ do-    a  <- m-    as <- interleaveSequence ms-    return (a:as)-{-# INLINE interleaveSequence #-}

@@ -8,295 +8,180 @@ --  module Control.Monad.MC.Sample (-    -- * Sampling from lists+    -- * Sampling+    -- ** Lists     sample,     sampleWithWeights,     sampleSubset,-    sampleSubset',     sampleSubsetWithWeights,-    sampleSubsetWithWeights',+    shuffle, -    -- * Sampling @Int@s+    -- ** Ints     sampleInt,     sampleIntWithWeights,     sampleIntSubset,-    sampleIntSubset',     sampleIntSubsetWithWeights,-    sampleIntSubsetWithWeights',--    -- * Shuffling-    shuffle,     shuffleInt,-    shuffleInt',     ) where -import Control.Monad-import Control.Monad.ST hiding (unsafeInterleaveST)-import Control.Monad.ST.Unsafe (unsafeInterleaveST)-import Control.Monad.MC.Base-import Control.Monad.MC.Repeat-import Control.Monad.MC.Walker+import Control.Monad( forM_, liftM )+import Control.Monad.Primitive( PrimMonad )+import Control.Monad.Trans.Class( lift ) import Data.List( foldl', sort )--import Data.Vector.Unboxed( MVector, Unbox ) import qualified Data.Vector as BV-import qualified Data.Vector.Mutable as BMV import qualified Data.Vector.Unboxed as V-import qualified Data.Vector.Generic.Mutable as MV+import qualified Data.Vector.Unboxed.Mutable as MV +import Control.Monad.MC.GSLBase+import Control.Monad.MC.Walker++ -- | @sample xs@ samples a value uniformly from the elements of @xs@.  The -- results are undefined if @length xs@ is zero.-sample :: (MonadMC m) => [a] -> m a+sample :: (PrimMonad m) => [a] -> MC m a sample xs = let     n = length xs     in sampleHelp n xs $ sampleInt n {-# INLINE sample #-} + -- | @sampleWithWeights wxs@ samples a value from the list with the given -- weight.-sampleWithWeights :: (MonadMC m) => [(Double, a)] -> m a+sampleWithWeights :: (PrimMonad m) => [(Double, a)] -> MC m a sampleWithWeights wxs = let     (ws,xs) = unzip wxs     n       = length xs     in sampleHelp n xs $ sampleIntWithWeights ws n {-# INLINE sampleWithWeights #-} ++sampleHelp :: (PrimMonad m) => Int -> [a] -> MC m Int -> MC m a+sampleHelp _n xs f = let+    arr = BV.fromList xs+    in liftM (BV.unsafeIndex arr) f+{-# INLINE sampleHelp #-}++ -- | @sampleSubset xs k@ samples a subset of size @k@ from @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.-sampleSubset :: (MonadMC m) => [a] -> Int -> m [a]+-- with the elements in the subset in the order that they were sampled.+sampleSubset :: (PrimMonad m) => [a] -> Int -> MC m [a] sampleSubset xs k = let     n = length xs-    in sampleListHelp n xs $ sampleIntSubset n k+    in sampleSubsetHelp n xs $ sampleIntSubset n k {-# INLINE sampleSubset #-} --- | Strict version of 'sampleSubset'.-sampleSubset' :: (MonadMC m) => [a] -> Int -> m [a]-sampleSubset' xs k = do-    s <- sampleSubset xs k-    length s `seq` return s-{-# INLINE sampleSubset' #-}  -- | Sample a subset of the elements with the given weights.  Return--- the elements of the subset lazily in the order they were sampled.-sampleSubsetWithWeights :: (MonadMC m) => [(Double,a)] -> Int -> m [a]+-- the elements of the subset in the order they were sampled.+sampleSubsetWithWeights :: (PrimMonad m) => [(Double,a)] -> Int -> MC m [a] sampleSubsetWithWeights wxs k = let     (ws,xs) = unzip wxs     n = length ws-    in sampleListHelp n xs $ sampleIntSubsetWithWeights ws n k+    in sampleSubsetHelp n xs $ sampleIntSubsetWithWeights ws n k {-# INLINE sampleSubsetWithWeights #-} --- | Strict version of 'sampleSubsetWithWeights'.-sampleSubsetWithWeights' :: (MonadMC m) => [(Double,a)] -> Int -> m [a]-sampleSubsetWithWeights' wxs k = do-    s <- sampleSubsetWithWeights wxs k-    length s `seq` return s-{-# INLINE sampleSubsetWithWeights' #-} -sampleHelp :: (Monad m) => Int -> [a] -> m Int -> m a-sampleHelp _n xs f = let-    arr = BV.fromList xs-    in liftM (BV.unsafeIndex arr) f-{-# INLINE sampleHelp #-}--sampleHelpU :: (Unbox a, Monad m) => Int -> [a] -> m Int -> m a-sampleHelpU _n xs f = let-    arr = V.fromList xs-    in liftM (V.unsafeIndex arr) f-{-# INLINE sampleHelpU #-}--{-# RULES "sampleHelp/Double" forall n xs f.-              sampleHelp n (xs :: [Double]) f = sampleHelpU n xs f #-}-{-# RULES "sampleHelp/Int" forall n xs f.-              sampleHelp n (xs :: [Int]) f = sampleHelpU n xs f #-}--sampleListHelp :: (Monad m) => Int -> [a] -> m [Int] -> m [a]-sampleListHelp _n xs f = let+sampleSubsetHelp :: (Monad m) => Int -> [a] -> m [Int] -> m [a]+sampleSubsetHelp _n xs f = let     arr = BV.fromList xs     in liftM (map $ BV.unsafeIndex arr) f-{-# INLINE sampleListHelp #-}+{-# INLINE sampleSubsetHelp #-} -sampleListHelpU :: (Unbox a, Monad m) => Int -> [a] -> m [Int] -> m [a]-sampleListHelpU _n xs f = let-    arr = V.fromList xs-    in liftM (map $ V.unsafeIndex arr) f -{-# RULES "sampleListHelp/Double" forall n xs f.-              sampleListHelp n (xs :: [Double]) f = sampleListHelpU n xs f #-}-{-# RULES "sampleListHelp/Int" forall n xs f.-              sampleListHelp n (xs :: [Int]) f = sampleListHelpU 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 :: (PrimMonad m) => Int -> MC 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 :: (PrimMonad m) => [Double] -> Int -> MC m Int sampleIntWithWeights ws n =     let qjs = computeTable n ws     in liftM (indexTable qjs) (uniform 0 1)-{-# INLINE sampleIntWithWeights #-} + -- | @sampleIntSubset n k@ 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]+-- were sampled.+sampleIntSubset :: (PrimMonad m) => Int -> Int -> MC m [Int] sampleIntSubset n k | k < 0     = fail "negative subset size"                     | k > n     = fail "subset size is too big"                     | otherwise = do-    us <- randomIndices n k-    return $ runST $ do-        ints <- MV.new n :: ST s (MVector s Int)-        sequence_ [ MV.unsafeWrite ints i i | i <- [0 .. n-1] ]-        sampleIntSubsetHelp ints us (n-1)+    xs <- lift $ (V.thaw . V.fromList) [ 0..n-1 ]+    go xs [] n k   where-    randomIndices n' k' | k' == 0   = return []-                        | otherwise = unsafeInterleaveMC $ do-        u  <- uniformInt n'-        us <- randomIndices (n'-1) (k'-1)-        return (u:us)+    go xs ys n' k' | k' == 0   = return $ reverse ys+                   | otherwise = do+        u <- uniformInt n'+        y <- lift $ do+                 i <- MV.unsafeRead xs u+                 j <- MV.unsafeRead xs (n'-1)+                 MV.unsafeWrite xs u j+                 return i+        go xs (y:ys) (n'-1) (k'-1) -    sampleIntSubsetHelp :: MVector m Int -> [Int] -> Int -> ST m [Int]-    sampleIntSubsetHelp _    []     _  = return []-    sampleIntSubsetHelp ints (u:us) n' = unsafeInterleaveST $ do-        i <- MV.unsafeRead ints u-        MV.unsafeWrite ints u =<< MV.unsafeRead ints n'-        is <- sampleIntSubsetHelp ints us (n'-1)-        return (i:is)-{-# INLINE sampleIntSubset #-} --- | Strict version of 'sampleIntSubset'.-sampleIntSubset' :: (MonadMC m) => Int -> Int -> m [Int]-sampleIntSubset' n k = do-    s <- sampleIntSubset n k-    length s `seq` return s-{-# INLINE sampleIntSubset' #-}- -- | @sampleIntSubsetWithWeights ws n k@ samplea size @k@ subset of -- @{ 0, ..., n-1 }@ with the given weights by sampling elements without--- replacement.  It returns the elements of the subset lazily in the order+-- replacement.  It returns the elements of the subset in the order -- they were sampled.-sampleIntSubsetWithWeights :: (MonadMC m) => [Double] -> Int -> Int -> m [Int]-sampleIntSubsetWithWeights ws n k = let-    w_sum0 = foldl' (+) 0 $ take n ws-    wjs = [ (w / w_sum0, j) | (w,j) <- reverse $ sort $ zip ws [ 0..n-1 ] ]+sampleIntSubsetWithWeights :: (PrimMonad m) => [Double] -> Int -> Int -> MC m [Int]+sampleIntSubsetWithWeights ws n k | k < 0     = fail "negative subset size"+                                  | k > n     = fail "subset size is too big"+                                  | otherwise = let+    wsum = foldl' (+) 0 $ take n ws+    wjs  = [ (w / wsum, j) | (w,j) <- reverse $ sort $ zip ws [ 0..n-1 ] ]     in do-        us <- replicateMC k $ uniform 0 1-        return $ runST $ do-            ints <- MV.new n :: ST s (MVector s (Double,Int))-            sequence_ [ MV.unsafeWrite ints i wj | (i,wj) <- zip [ 0.. ] wjs ]-            go ints n 1 us+        xs <- lift $ (V.thaw . V.fromList) wjs+        go xs wsum [] n k   where-    go :: MVector m (Double, Int) -> Int -> Double -> [Double] -> ST m [Int]-    go ints n' w_sum us | null us   = return []-                        | otherwise = let-        target = head us * w_sum-        in unsafeInterleaveST $ do-            (i,(w,j)) <- findTarget ints n' target 0 0-            shiftDown ints (i+1) (n'-1)-            let w_sum' = w_sum - w-                n''    = n' - 1-                us'    = tail us-            js <- go ints n'' w_sum' us'-            return $ j:js+    go xs wsum' ys n' k' | k' == 0   = return $ reverse ys+                         | otherwise = do+        target <- uniform 0 wsum'+        (w,y) <- lift $ do+                     (i,wj) <- findTarget xs n' target 0 0+                     shiftDown xs (i+1) (n'-1)+                     return wj+        let wsum'' = wsum' - w+            ys'    = y:ys+            n''    = n' - 1+            k''    = k' - 1+        go xs wsum'' ys' n'' k'' -    findTarget :: MVector m (Double, Int)-                    -> Int -> Double -> Int -> Double -> ST m (Int, (Double, Int))-    findTarget ints n' target i acc+    findTarget xs n' target i acc         | i == n' - 1 = do-            wj <- MV.unsafeRead ints i+            wj <- MV.unsafeRead xs i             return (i,wj)         | otherwise = do-            (w,j) <- MV.unsafeRead ints i+            (w,j) <- MV.unsafeRead xs i             let acc' = acc + w             if target <= acc'                 then return (i,(w,j))-                else findTarget ints n' target (i+1) acc'+                else findTarget xs n' target (i+1) acc' -    shiftDown :: MVector m (Double, Int) -> Int -> Int -> ST m ()-    shiftDown ints from to =+    shiftDown xs from to =         forM_ [ from..to ] $ \i -> do-            wj <- MV.unsafeRead ints i-            MV.unsafeWrite ints (i-1) wj+            wj <- MV.unsafeRead xs i+            MV.unsafeWrite xs (i-1) wj -{-# INLINE sampleIntSubsetWithWeights #-} --- | Strict version of 'sampleIntSubsetWithWeights'.-sampleIntSubsetWithWeights' :: (MonadMC m) => [Double] -> Int -> Int -> m [Int]-sampleIntSubsetWithWeights' ws n k = do-    s <- sampleIntSubsetWithWeights ws n k-    length s `seq` return s-{-# INLINE sampleIntSubsetWithWeights' #-}- -- | @shuffle xs@ randomly permutes the list @xs@ and returns -- the result.  All permutations of the elements of @xs@ are equally -- likely.-shuffle :: (MonadMC m) => [a] -> m [a]+shuffle :: (PrimMonad m) => [a] -> MC m [a] shuffle xs = let-    n = length xs-    in shuffleInt n >>= \swaps -> (return . BV.toList) $ BV.create $ do-           marr <- MV.new n-           zipWithM_ (MV.unsafeWrite marr) [0 .. n-1] xs-           mapM_ (swap marr) swaps-           return marr-  where-    swap :: BMV.MVector s a -> (Int,Int) -> ST s ()-    swap marr (i,j) | i == j    = return ()-                    | otherwise = do-        x <- MV.unsafeRead marr i-        y <- MV.unsafeRead marr j-        MV.unsafeWrite marr i y-        MV.unsafeWrite marr j x-{-# INLINE shuffle #-}--shuffleU :: (Unbox a, MonadMC m) => [a] -> m [a]-shuffleU xs = let-    n = length xs-    in shuffleInt n >>= \swaps -> (return . V.toList) $ V.create $ do-           marr <- MV.new n-           zipWithM_ (MV.unsafeWrite marr) [0 .. n-1] xs-           mapM_ (swap marr) swaps-           return marr-  where-    swap :: (Unbox a) => MVector s a -> (Int,Int) -> ST s ()-    swap marr (i,j) | i == j    = return ()-                    | otherwise = do-        x <- MV.unsafeRead marr i-        y <- MV.unsafeRead marr j-        MV.unsafeWrite marr i y-        MV.unsafeWrite marr j x-{-# INLINE shuffleU #-}--{-# RULES "shuffle/Double" forall xs.-              shuffle (xs :: [Double]) = shuffleU xs #-}-{-# RULES "shuffle/Int" forall xs.-              shuffle (xs :: [Int]) = shuffleU xs #-}+    n   = length xs+    mis = liftM BV.fromList $ shuffleInt n+    in liftM (BV.toList . BV.unsafeBackpermute (BV.fromList xs)) mis  --- | @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 #-}---- | Strict version of 'shuffleInt'.-shuffleInt' :: (MonadMC m) => Int -> m [(Int,Int)]-shuffleInt' n = do-    ss <- shuffleInt n-    length ss `seq` return ss-{-# INLINE shuffleInt' #-}+-- | @shuffleInt n@ randomly permutes the elements of the list @[ 0..n-1 ]@.+shuffleInt :: (PrimMonad m) => Int -> MC m [Int]+shuffleInt n = sampleIntSubset n n

@@ -1,15 +0,0 @@--------------------------------------------------------------------------------- |--- Module     : Data.Summary--- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>--- License    : BSD3--- Maintainer : Patrick Perry <patperry@gmail.com>--- Stability  : experimental------ Summary Statistics----module Data.Summary (-    module Data.Summary.Double-    ) where--import Data.Summary.Double

@@ -1,3 +1,4 @@+{-# LANGUAGE DeriveDataTypeable #-} ----------------------------------------------------------------------------- -- | -- Module     : Data.Summary.Bool@@ -10,26 +11,48 @@ --  module Data.Summary.Bool (-    -- * The @Summary@ data type+    -- * Summary type     Summary,-    summary,-    update, -    -- * @Summary@ properties-    sampleSize,-    count,-    sampleMean,-    sampleSE,-    sampleCI,+    -- * Properties+    -- ** Sum+    size,+    sum,+    -- ** Mean+    mean,+    meanSE,+    meanCI, +    -- * Construction+    empty,+    singleton,++    -- * Insertion+    insert,+    insertWith,++    -- * Combination+    union,+    unions,++    -- * Conversion+    -- ** Lists+    fromList,+    fromListWith,++    -- ** Statistics+    toStats,+    fromStats,+     ) where -import Control.DeepSeq+import Prelude hiding (sum) import Data.List( foldl' )-import Data.Monoid-import Text.Printf+import Data.Monoid( Monoid(..) )+import Data.Data( Data, Typeable )+import Text.Printf( printf ) -import Data.Summary.Utils+import Data.Summary.Utils( interval )   -- | A type for storing summary statistics for a data set of@@ -37,63 +60,86 @@ -- of 'True' events and gives estimates for the success -- probability.  'True' is interpreted as a one, and 'False' -- is interpreted as a zero.-data Summary = S {-# UNPACK #-} !Int  -- sample size-                 {-# UNPACK #-} !Int  -- number of successes+data Summary = S {-# UNPACK #-} !Int -- number of observations+                 {-# UNPACK #-} !Int -- number of True values+    deriving(Eq, Data, Typeable)  instance Show Summary where-    show s@(S n c) =-        printf "    sample size: %d" n-        ++ printf "\n      successes: %d" c-        ++ printf "\n     proportion: %g" (sampleMean s)-        ++ printf "\n             SE: %g" (sampleSE s)+    show s@(S n x) =+             printf "    sample size: %d" n+        ++ printf "\n      successes: %d" x+        ++ printf "\n     proportion: %g" (mean s)+        ++ printf "\n             SE: %g" (meanSE s)         ++ printf "\n         99%% CI: (%g, %g)" c1 c2-      where (c1,c2) = sampleCI 0.99 s+      where (c1,c2) = meanCI 0.99 s  instance Monoid Summary where     mempty = empty     mappend = union -instance NFData Summary+-- | Number of observations.+size :: Summary -> Int+size (S n _) = n +-- | Number of 'True' values.+sum :: Summary -> Int+sum (S _ x) = x --- | Get a summary of a list of values.-summary :: [Bool] -> Summary-summary = foldl' update empty+-- | Proportion of 'True' values.+mean :: Summary -> Double+mean (S n x) = fromIntegral x / fromIntegral n +-- | Standard error for the mean (proportion of 'True' values).+meanSE :: Summary -> Double+meanSE s = sqrt (p*(1-p) / n)+  where p = mean s+        n = fromIntegral $ size s++-- | Central Limit Theorem based confidence interval for the+-- population mean (proportion) at the specified coverage level.  The+-- level must be in the range @(0,1)@.+meanCI :: Double -> Summary -> (Double,Double)+meanCI level s = interval level (mean s) (meanSE s)+ -- | Get an empty summary. empty :: Summary empty = S 0 0 --- | Take the union of two summaries.-union :: Summary -> Summary -> Summary-union (S na ca) (S nb cb) = S (na + nb) (ca + cb)+-- | Summarize a single value.+singleton :: Bool -> Summary+singleton x = S 1 (if x then 1 else 0)  -- | Update the summary with a data point.-update :: Summary -> Bool -> Summary-update (S n c) i =+insert :: Bool -> Summary -> Summary+insert y (S n x) =     let n' = n+1-        c' = if i then c+1 else c-    in S n' c'+        x' = if y then x+1 else x+    in S n' x' --- | Get the sample size.-sampleSize :: Summary -> Int-sampleSize (S n _) = n+-- | Apply a function and update the summary with the result.+insertWith :: (a -> Bool) -> a -> Summary -> Summary+insertWith f a = insert (f a) --- | Get the number of 'True' values-count :: Summary -> Int-count (S _ c) = c+-- | Take the union of two summaries.+union :: Summary -> Summary -> Summary+union (S na xa) (S nb xb) = S (na + nb) (xa + xb) --- | Get the proportion of 'True' events.-sampleMean :: Summary -> Double-sampleMean (S n c) = fromIntegral c / fromIntegral n+-- | Take the union of a list of summaries.+unions :: [Summary] -> Summary+unions = foldl' union empty --- | Get the standard error for the sample proportion.-sampleSE :: Summary -> Double-sampleSE s = sqrt (p*(1-p) / n)-  where p = sampleMean s-        n = fromIntegral $ sampleSize s+-- | Get a summary of a list of values.+fromList :: [Bool] -> Summary+fromList = foldl' (flip insert) empty --- | 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 = interval level (sampleMean s) (sampleSE s)+-- | Map a function over a list of values and summarize the results.+fromListWith :: (a -> Bool) -> [a] -> Summary+fromListWith f = fromList . map f++-- | Convert to (size,sum).+toStats :: Summary -> (Int,Int)+toStats (S n x) = (n,x)++-- | Convert from (size,sum).  No validation is performed.+fromStats :: Int -> Int -> Summary+fromStats = S

@@ -1,3 +1,4 @@+{-# LANGUAGE DeriveDataTypeable #-} ----------------------------------------------------------------------------- -- | -- Module     : Data.Summary.Double@@ -10,29 +11,57 @@ --  module Data.Summary.Double (-    -- * The @Summary@ data type+    -- * Summary type     Summary,-    summary,-    update, -    -- * @Summary@ properties-    sampleSize,-    sampleMin,-    sampleMax,-    sampleMean,-    sampleSE,-    sampleVar,-    sampleSD,-    sampleCI,+    -- * Properties+    -- ** Sum+    size,+    sum,+    sumSquaredErrors,+    -- ** Mean+    mean,+    meanSE,+    meanCI,+    -- ** Scale+    stddev,+    variance,+    -- ** Extremes+    min,+    max, +    -- * Construction+    empty,+    singleton,++    -- * Insertion+    insert,+    insertWith,++    -- * Combination+    union,+    unions,++    -- * Conversion+    -- ** Lists+    fromList,+    fromListWith,++    -- ** Statistics+    toStats,+    fromStats,+     ) where -import Control.DeepSeq+import Prelude hiding (min, max, sum)+import qualified Prelude as P++import Data.Data( Data, Typeable ) import Data.List( foldl' ) import Data.Monoid import Text.Printf -import Data.Summary.Utils+import Data.Summary.Utils( interval )   -- | A type for storing summary statistics for a data set including@@ -42,38 +71,78 @@                  {-# UNPACK #-} !Double  -- sum of squares                  {-# UNPACK #-} !Double  -- sample min                  {-# UNPACK #-} !Double  -- sample max+    deriving(Eq, Data, Typeable)  instance Show Summary where     show s@(S n mu _ l h) =-        printf "    sample size: %d" n+             printf "    sample size: %d" n         ++ printf "\n            min: %g" l         ++ printf "\n            max: %g" h         ++ printf "\n           mean: %g" mu-        ++ printf "\n             SE: %g" (sampleSE s)+        ++ printf "\n             SE: %g" (meanSE s)         ++ printf "\n         99%% CI: (%g, %g)" c1 c2-      where (c1,c2) = sampleCI 0.99 s+      where (c1,c2) = meanCI 0.99 s  instance Monoid Summary where     mempty = empty     mappend = union -instance NFData Summary+-- | Number of observations.+size :: Summary -> Int+size (S n _ _ _ _) = n +-- | Sum of values.+sum :: Summary -> Double+sum s = (fromIntegral $ size s) * (mean s) --- | Get a summary of a list of values.-summary :: [Double] -> Summary-summary = foldl' update empty+-- | Mean value.+mean :: Summary -> Double+mean (S _ m _ _ _) = m --- | Get an empty summary.+-- | Standard error of the mean.+meanSE :: Summary -> Double+meanSE s = sqrt (variance s / fromIntegral (size s))++-- | Get a Central Limit Theorem based confidence interval for the+-- population mean with the specified coverage level.  The level must+-- be in the range @(0,1)@.+meanCI :: Double -> Summary -> (Double,Double)+meanCI level s = interval level (mean s) (meanSE s)++-- | Sample standard deviation.+stddev :: Summary -> Double+stddev s = sqrt (variance s)++-- | Sample variance.+variance :: Summary -> Double+variance s = (sumSquaredErrors s) / fromIntegral (size s - 1)++-- | Sum of squared errors @(x[i] - mean)^2@.+sumSquaredErrors :: Summary -> Double+sumSquaredErrors (S _ _ s _ _) = s++-- | Minimum value.+min :: Summary -> Double+min (S _ _ _ l _) = l++-- | Maximum value.+max :: Summary -> Double+max (S _ _ _ _ h) = h++-- | An empty summary. empty :: Summary empty = S 0 0 0 (1/0) (-1/0) +-- | Summarize a single value.+singleton :: Double -> Summary+singleton x = S 1 x 0 x x+ -- | 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+-- 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 =+insert :: Double -> Summary -> Summary+insert x (S n m s l h) =     let n'    = n+1         delta = x - m         m'    = m + delta / fromIntegral n'@@ -82,10 +151,14 @@         h'    = if x > h then x else h     in S n' m' s' l' h' +-- | Apply a function and update the summary with the result.+insertWith :: (a -> Double) -> a -> Summary -> Summary+insertWith f a = insert (f a)+ -- | Take the union of two summaries.--- Use the updating rules from Chan et al. "Updating Formulae and a Pairwise---   Algorithm for Computing Sample Variances," available at--- ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/79/773/CS-TR-79-773.pdf+-- Use the updating rules from Chan et al. \"Updating Formulae and a Pairwise+--   Algorithm for Computing Sample Variances,\" available at+--   <http://infolab.stanford.edu/pub/cstr/reports/cs/tr/79/773/CS-TR-79-773.pdf>. union :: Summary -> Summary -> Summary union (S na ma sa la ha) (S nb mb sb lb hb) =     let delta = mb - ma@@ -97,39 +170,28 @@            | otherwise = ma + weightedDelta         s  | n == 0    = 0            | otherwise = sa + sb + delta*na'*weightedDelta-        l  = min la lb-        h  = max ha hb+        l  = P.min la lb+        h  = P.max ha hb     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)+-- | Take the union of a list of summaries.+unions :: [Summary] -> Summary+unions = foldl' union empty --- | Get the sample standard deviation.-sampleSD :: Summary -> Double-sampleSD s = sqrt (sampleVar s)+-- | Get a summary of a list of values.+fromList :: [Double] -> Summary+fromList = foldl' (flip insert) empty --- | Get the sample standard error.-sampleSE :: Summary -> Double-sampleSE s = sqrt (sampleVar s / fromIntegral (sampleSize s))+-- | Map a function over a list of values and summarize the results.+fromListWith :: (a -> Double) -> [a] -> Summary+fromListWith f = fromList . map f --- | 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 = interval level (sampleMean s) (sampleSE s)+-- | Convert to (size, mean, sumSquaredErrors, min, max).+toStats :: Summary -> (Int,Double,Double,Double,Double)+toStats (S n m s l u) = (n,m,s,l,u) --- | Get the minimum of the sample.-sampleMin :: Summary -> Double-sampleMin (S _ _ _ l _) = l+-- | Convert from (size, mean, sumSquaredErrors, min, max).+-- No validation is performed.+fromStats :: Int -> Double -> Double -> Double -> Double -> Summary+fromStats = S --- | Get the maximum of the sample.-sampleMax :: Summary -> Double-sampleMax (S _ _ _ _ h) = h

@@ -16,12 +16,13 @@  import GSL.Random.Dist( ugaussianPInv ) --- | Get a Central Limit Theorem-based confidence interval for the--- population mean with the specified coverage level.  The level must--- be in the range @(0,1)@.-interval :: Double -- ^ the confidence level-         -> Double -- ^ the sample mean-         -> Double -- ^ the sample standard error++-- | Get a Central Limit Theorem based confidence interval for a+-- population parameter with the specified coverage level.  The level+-- must be in the range @(0,1)@.+interval :: Double -- ^ Confidence level.+         -> Double -- ^ Estimate.+         -> Double -- ^ Standard error.          -> (Double,Double) interval level xbar se | not (level > 0 && level < 1) =                              error "level must be between 0 and 1"@@ -31,6 +32,7 @@         delta = z*se     in (xbar-delta, xbar+delta) --- | Tests if the value is in the open interval (a,b)++-- | Tests if the value is in the open interval (a,b). inInterval :: Double -> (Double,Double) -> Bool x `inInterval` (a,b) = x > a && x < b

@@ -1,5 +1,5 @@ name:           monte-carlo-version:        0.5+version:        0.6 license:        BSD3 license-file:   LICENSE author:         Patrick Perry@@ -7,47 +7,45 @@ homepage:       http://github.com/patperry/hs-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.+description:    A monad and transformer for performing Monte Carlo+                calculations. This monad carries and provides access to+                a pseudo-random number generator. Internally, the monad+                mutates rather than copies the random gnerator state.  By+                avoiding copies, it can deliver faster performance than+                many pure random number implementations.  The package is+                built around the facilities provided by the GNU Scientific+                Library (GSL). build-type:     Simple stability:      experimental cabal-version:  >= 1.8 extra-source-files: NEWS examples/Binomial.hs examples/Pi.lhs-                    examples/Poker.hs  examples/Queue.hs tests/Main.hs+                    examples/Poker.hs examples/Queue.hs tests/Main.hs                     tests/Makefile  library     exposed-modules:             Control.Monad.MC-            Control.Monad.MC.Class-            Control.Monad.MC.GSL-            Data.Summary             Data.Summary.Bool             Data.Summary.Double-            Data.Summary.Utils      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.Walker+            Data.Summary.Utils      extensions:-            FlexibleInstances,-            MultiParamTypeClasses,-            TypeFamilies,-            UndecidableInstances+            DeriveDataTypeable+            RankNTypes+            TypeFamilies -    build-depends:  base       >= 4     && < 5,-                    gsl-random >= 0.4.3 && < 0.5,-                    mtl        >= 1.1   && < 3.0,-                    vector     >= 0.6   && < 0.11,-                    deepseq    >= 1.0   && < 2.0+    build-depends:  base         >= 4     && < 5,+                    gsl-random   >= 0.5   && < 0.6,+                    primitive    >= 0.5   && < 0.6,+                    transformers >= 0.3   && < 0.5,+                    vector       >= 0.6   && < 0.11      hs-source-dirs: lib     ghc-options:    -Wall@@ -66,13 +64,14 @@   ghc-options:       -Wall -Werror   build-depends:-      base       >= 4     && < 5-    , gsl-random >= 0.4.3 && < 0.5-    , mtl        >= 1.1   && < 3.0-    , vector     >= 0.6   && < 0.11-    , ieee754 >= 0.7 && < 0.8-    , random >=1.0 && < 1.1-    , QuickCheck >= 2.4.0.1 && < 2.6-    , test-framework-quickcheck2 >= 0.2 && < 0.4-    , test-framework >= 0.4 && < 0.9-    , deepseq >= 1.0 && < 2.0+      base                       >= 4       && < 5,+      gsl-random                 >= 0.5     && < 0.6,+      primitive                  >= 0.5     && < 0.6,+      transformers               >= 0.3     && < 0.5,+      vector                     >= 0.6     && < 0.11,+      ieee754                    >= 0.7     && < 0.8,+      random                     >= 1.0     && < 1.1,+      QuickCheck                 >= 2.4.0.1 && < 2.8,+      test-framework             >= 0.4     && < 0.9,+      test-framework-quickcheck2 >= 0.2     && < 0.4+

@@ -4,7 +4,8 @@ import Control.Monad import Data.AEq import Data.Monoid-import Data.Summary+import Data.Summary.Double( Summary )+import qualified Data.Summary.Double as S import Test.QuickCheck import Test.Framework import Test.Framework.Providers.QuickCheck2@@ -40,17 +41,17 @@  prop_monoid_update_equiv :: [Double] -> [Double] -> Bool prop_monoid_update_equiv xs ys =-    approxEqualS (summary $ xs <> ys)-                 (summary xs <> summary ys)+    approxEqualS (S.fromList $ xs <> ys)+                 (S.fromList xs <> S.fromList ys)  prop_monoid_assoc :: [Double] -> [Double] -> [Double] -> Bool prop_monoid_assoc xs ys zs =-    let (sxs, sys, szs) = (summary xs, summary ys, summary zs)+    let (sxs, sys, szs) = (S.fromList xs, S.fromList ys, S.fromList zs)      in ((sxs <> sys) <> szs) `approxEqualS` (sxs <> (sys <> szs))  prop_monoid_commute :: [Double] -> [Double] -> Bool prop_monoid_commute xs ys =-    let (sxs, sys) = (summary xs, summary ys)+    let (sxs, sys) = (S.fromList xs, S.fromList ys)      in (sxs <> sys) `approxEqualS` (sys <> sxs)  tests_monoid :: Test@@ -70,8 +71,8 @@  approxEqualS :: Summary -> Summary -> Bool approxEqualS a b =-    sampleSize a == sampleSize b &&-      all eq [ sampleMin, sampleMax, sampleMean, sampleVar ]+    S.size a == S.size b &&+      all eq [ S.min, S.max, S.mean, S.variance ]     where         eq f = f a ~== f b