I’d like to know if I’m doing this the most performant way. I’m pretty new to Haskell so I’m not very familiar with the ecosystem or performance considerations yet.
I need to generate some random data representing how a marketing campaign is doing in terms of how successfully the content is being delivered. For example, there are push notifications, email, etc. For each of these channels, a campaign could succeed or fail to be delivered. So for any given campaign, there will b N number emails sent, and some subset of them may fail. However, the randomness of success/failure is not evenly distributed. Each campaign has a predefined probability of experiencing failure, so that probability is factored into the random outcomes. The way I’m currently generating those outcomes is like this:
module Main where
import qualified Control.Monad.Random as CMR
import Data.Function ( (&) )
genFailures :: Int -> IO Int
genFailures numSent =
-- True is Success; False is Failure
CMR.weighted [ (True, 8 / 10), (False, 2 / 10) ]
& CMR.replicateM numSent
& fmap (length . filter (== False))
main :: IO ()
main =
genFailures 1000 >>= print
When numSent is large, this is very slow and I’m wondering if there is perhaps a more performant way of achieving what I want. Any pointers would be much appreciated!
A more immediate way to improve your code is to keep track of the number of failures while doing the random number generation:
{-# LANGUAGE BangPatterns #-}
module Main where
import qualified Control.Monad.Random as CMR
genFailures :: Int -> IO Int
genFailures n = genFailures' 0 n
genFailures' :: Int -> Int -> IO Int
genFailures' acc 0 = return acc
genFailures' !acc !n = do
x <- CMR.weighted [ (True, 8 / 10), (False, 2 / 10) ]
genFailures' (if x then acc else acc + 1) (n - 1)
main :: IO ()
main = genFailures 1000000 >>= print
Just tried using mwc-probability (System.Random.MWC.Probability) for binomial sampling, but it appears to be even slower than my initial solution with MonadRandom (Control.Monad.Random). I’m testing both again 1_000_000 as input.
$ time stack exec monad-random-gen-failures-exe
> 1.83 real 2.28 user 0.31 sys
$ time stack exec mwc-prob-gen-failures-exe
> 1.97 real 5.95 user 1.38 sys
Looks like this:
module Main where
import qualified System.Random.MWC.Probability as Prob
binomialFailures :: Prob.GenIO -> Int -> IO Int
binomialFailures gen numSent =
Prob.sample (Prob.binomial numSent (4 / 5)) gen
main :: IO ()
main = do
gen <- Prob.createSystemRandom
binomialFailures gen 1000000
>>= print
I figured it might be faster since MonadRandom is based on random (System.Random) and mwc-probability is based on mwc-random.
@jaror thanks for your input. For some reason, when I tried your recursive version, it was slower than my original. My types are little different but I don’t see why it would affect the performance. Looks like this:
{-# LANGUAGE BangPatterns #-}
module Main where
import qualified Control.Monad.Random as CMR
genFailures :: Int -> CMR.Rand CMR.StdGen Int
genFailures =
genFailures' 0
genFailures' :: Int -> Int -> CMR.Rand CMR.StdGen Int
genFailures' acc 0 =
return acc
genFailures' !acc !n =
do
x <- CMR.weighted [ (True, 8 / 10), (False, 2 / 10) ]
genFailures' (if x then acc else acc + 1) (n - 1)
main :: IO ()
main = CMR.evalRandIO (genFailures 1000000) >>= print
Here are the times I’m seeing:
$ time stack exec original-exe
> 1.81 real 2.35 user 0.34 sys
$ time stack exec recursive-exe
> 3.16 real 3.88 user 0.49 sys
Yes, I noticed after writing it that my version was not that much more efficient than yours. I think that the Rand StdGen Int type and the evalRandIO are the culprits of the slowdown. But it is probably not efficient enough anyway.
mwc-random is also my go-to package if I need fast-ish random numbers. If we just use a naive tail-recursive version, we get:
module Main where
import Control.Monad.Primitive (PrimMonad (..))
import qualified System.Random.MWC as MWC
genFailures :: PrimMonad m => MWC.Gen (PrimState m) -> Int -> m Int
genFailures gen = go 0
where
go !acc n
| n <= 0 = return acc
| otherwise = do
x <- MWC.uniformR (1, 10 :: Int) gen
go (if x <= 8 then acc else acc + 1) (n - 1)
main :: IO ()
main = do
gen <- MWC.createSystemRandom
print =<< genFailures gen 1000000
Which gives us a decent speedup:
> time ./original
200431
real 0m2.254s
user 0m2.127s
sys 0m0.119s
> time ./mwc
200091
real 0m0.039s
user 0m0.034s
sys 0m0.005s
However, the code is very imperative! Two common problems with the solutions posted so far are:
weighted seems slow and doesn’t get inlined properly.
replicateM in IO is not lazy enough – it will generate the entire list before starting to count Falses.
We can write a more functional version using System.Random easily:
{-# LANGUAGE BangPatterns #-}
module Main where
import System.Random
genFailures :: Int -> StdGen -> Int
genFailures n = length . filter (> 8) . take n . randomRs (1, 10 :: Int)
main :: IO ()
main = do
gen <- getStdGen
print $ genFailures 1000000 gen
But it’s a bit slower again, due to the use of StdGen over mwc-random:
> time ./functional
199713
real 0m0.198s
user 0m0.191s
sys 0m0.007s
If you want to avoid the imperative style but still want the speed, you can probably use (unsafe) interleaving to generate a lazy list of random numbers using mwc-random.
Thanks for your help @jaspervdj! Does your last example using System.Random account for my requirement that each campaign has a predefined probability of success that be used as a weight when generating the random outcomes? It doesn’t seem like it does, but I may just be missing something.
I’ve tried several different options and they’ve all been unacceptably slow for my purposes. I ended up using a hack to approximate a random sampling from a binomial distribution. The basic idea is to randomly add or subtract a random number that is between 1 and 3 standard deviations from the average using the 68-95-99.7 rule.
module Main where
import Data.Function ( (&) )
import qualified Control.Monad.Random as CMR
import qualified Statistics.Distribution as Stats
import qualified Statistics.Distribution.Binomial as Stats
main :: IO ()
main = CMR.evalRandIO (randErrCountFromTotal 1000000 0.5) >>= print
randErrCountFromTotal :: (CMR.RandomGen g) => Int -> Double -> CMR.Rand g Int
randErrCountFromTotal numEvts successRate = do
-- add or subtract 1-3 std deviations from the average
-- using the 68–95–99.7 rule as a heuristic
numStdDeviations <- CMR.fromList [ (1.0, 0.6827), (2.0, 0.2718), (3.0, 0.0455) ]
let stdDeviation = Stats.stdDev $ Stats.binomial numEvts successRate
deviationMin = round $ (numStdDeviations * stdDeviation) - stdDeviation
deviationMax = round $ numStdDeviations * stdDeviation
avgSuccessCount = round $ realToFrac numEvts * successRate
CMR.fromList [ ((+), 1 / 2), ((-), 1 / 2) ]
<*> pure avgSuccessCount
<*> CMR.getRandomR (deviationMin, deviationMax)
& fmap (\successCount -> numEvts - successCount)
Not being a math person myself, I’m not sure if this is all that sound of an approach, but it seems to approximate the effect I want pretty well and it’s much faster than the alternatives I tried.
$ time stack exec hacky-exe
> 0.27 real 0.21 user 0.03 sys
I think it really depends on what you will be using this for. It is pretty easy (i.e. not a lot of samples are required) to detect that your distribution is not actually a binomial distribution. But the samples are generally in the right place. It is like a really pixelated image of a nice continuous curve. It is like most things in statistics: it is not perfect, but it might be good enough for your purpose.
In this case I really think that it would be worthwhile to translate the algorithm that R uses to Haskell or maybe even just bind the C file. I have modified the rbinom.c file for easier binding here: https://gist.github.com/noughtmare/ff3705c0906506911266b45d03c11d20.
I have not tested it extensively, but it seems to work.
Now you can simulate running a marketing campaign on the whole population of earth multiple times over in 0 seconds, using main = print =<< runRandIO (rbinom 10000000000 0.2):
$ cabal new-run rbinom.hs
Resolving dependencies...
Build profile: -w ghc-8.6.5 -O1
In order, the following will be built (use -v for more details):
- fake-package-0 (exe:script) (first run)
Configuring executable 'script' for fake-package-0..
Preprocessing executable 'script' for fake-package-0..
Building executable 'script' for fake-package-0..
[1 of 1] Compiling Main ( Main.hs, /tmp/cabal-repl.-2130/dist-newstyle/build/x86_64-linux/ghc-8.6.5/fake-package-0/x/script/build/script/script-tmp/Main.o )
Linking /tmp/cabal-repl.-2130/dist-newstyle/build/x86_64-linux/ghc-8.6.5/fake-package-0/x/script/build/script/script ...
1999924842
93,344 bytes allocated in the heap
4,376 bytes copied during GC
53,688 bytes maximum residency (1 sample(s))
28,232 bytes maximum slop
0 MB total memory in use (0 MB lost due to fragmentation)
Tot time (elapsed) Avg pause Max pause
Gen 0 0 colls, 0 par 0.000s 0.000s 0.0000s 0.0000s
Gen 1 1 colls, 0 par 0.000s 0.000s 0.0001s 0.0001s
INIT time 0.000s ( 0.000s elapsed)
MUT time 0.000s ( 0.000s elapsed)
GC time 0.000s ( 0.000s elapsed)
EXIT time 0.000s ( 0.000s elapsed)
Total time 0.000s ( 0.000s elapsed)
%GC time 0.0% (0.0% elapsed)
Alloc rate 462,099,009 bytes per MUT second
Productivity 50.6% of total user, 51.6% of total elapsed
EDIT: I have tested it some more and it seems that it deviates to smaller values.