Skip to content

Instantly share code, notes, and snippets.

# dimitri-xyz/RejectionSampling.hs Created Apr 20, 2017

A non-biased way to generate random integers in a given range
 import Control.Monad.Random.Class import System.Random import Data.Word maxIterations :: Int maxIterations = 100 acceptReject :: (Random a, MonadRandom m) => (a -> Bool) -> (a -> b) -> m b acceptReject = acceptReject' (Just maxIterations) -- First argument may specify max number of iterations acceptReject' :: (Random a, MonadRandom m) => Maybe Int -> (a -> Bool) -> (a -> b) -> m b acceptReject' (Just 0) _ _ = error "acceptReject': Too many failed iterations, check your source of randomness." acceptReject' iter accept convert = do sample <- getRandom if accept sample then return (convert sample) else acceptReject' (subtract 1 <\$> iter) accept convert -- maxBound is 1 less than what we want, but that's ok, -- it will just make the range (maxBound+1)/2 be innefficiently mapped, but it will only be -- as slow as (maxBound+1)/2 + 1 already was (i.e. about twice as many iterations needed) inUniformRange :: Word64 -> Word64 -> Bool inUniformRange i x = x < (maxBound `div` i) * i constrainToRange :: Word64 -> Word64 -> Word64 constrainToRange range sample = sample `mod` range getRandLessThan :: MonadRandom m => Word64 -> m Word64 getRandLessThan range = acceptReject (inUniformRange range) (constrainToRange range) -- This may fail if: -- h-l+1 > (maxBound :: Int) -- OR h-l+1 > (maxBound :: Word64) getRandInRange :: MonadRandom m => Int -> Int -> m Int getRandInRange l h | l > h = error "getRandInRange: lower bound must be smaller than or equal to upper bound" getRandInRange l h = do r <- getRandLessThan (fromIntegral \$ h - l + 1) return \$ l + fromIntegral r
to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.