diff --git a/System/Random/MWC/Distributions.hs b/System/Random/MWC/Distributions.hs --- a/System/Random/MWC/Distributions.hs +++ b/System/Random/MWC/Distributions.hs @@ -29,6 +29,7 @@ , geometric1 , bernoulli , binomial + , poisson -- ** Multivariate , dirichlet -- * Permutations @@ -45,7 +46,7 @@ import Data.Foldable (foldl') import Data.Traversable (mapM) import Data.Word (Word32) -import System.Random.Stateful (StatefulGen(..),Uniform(..),UniformRange(..),uniformDoublePositive01M) +import System.Random.Stateful (StatefulGen(..),Uniform(..),UniformRange(..),uniformDoublePositive01M, uniformDouble01M) import qualified Data.Vector.Unboxed as I import qualified Data.Vector.Generic as G import qualified Data.Vector.Generic.Mutable as M @@ -353,6 +354,8 @@ -- \[ -- f(k;n,p) = \Pr(X = k) = \binom n k p^k(1-p)^{n-k} -- \] +-- +-- @since 0.15.1.0 binomial :: forall g m . StatefulGen g m => Int -- ^ Number of trials, must be positive. -> Double -- ^ Probability of success \(p \in [0,1]\) @@ -484,6 +487,81 @@ r' = r * ((a / fromIntegral x') - s) +-- | Random variate generate for Poisson distribution. +-- +-- If parameter λ is within 10 σ or greater than @maxBound :: Int@ +-- error is raised since result may not be representable as @Int@ +-- +-- @since 0.15.3.0 +poisson + :: StatefulGen g m + => Double -- ^ Rate parameter, also known as \( \lambda \) + -> g -- ^ Generator + -> m Int +{-# INLINE poisson #-} +poisson lambda gen + | lambda > maxPoissonLam + = pkgError "poisson" + $ "lambda is too large: "++show lambda++" Result won't fit into Int" + | lambda > 10 = poissonAtkinson lambda gen + | lambda >= 0 = poissonInterArrival lambda gen + | otherwise + = pkgError "poisson" "Lambda parameter must be greater than zero" + +maxPoissonLam :: Double +maxPoissonLam = m - 10 * sqrt m where + m = fromIntegral (maxBound :: Int) + + +-- This uses the fact that if N(t) is a Poisson process +-- with rate lambda, then the counting process N(t) can +-- be represented as interarrival times X[1], X[2],... with +-- X[i] ~ Exp(lambda). +poissonInterArrival :: StatefulGen g m + => Double -- ^ Rate parameter, also known as lambda + -> g -- ^ Generator + -> m Int +{-# INLINE poissonInterArrival #-} +poissonInterArrival lambda gen = do + loop 0 1.0 + where + loop !k !p = do + p' <- (*p) <$> uniformDouble01M gen + if p' <= bigL then return $! k else loop (k+1) p' + bigL = exp (negate lambda) + +-- Attributed to Atkinson, via Casella. Uses a rejection +-- method that uses logistic distribution as the envelope +-- distribution. +poissonAtkinson :: forall g m . StatefulGen g m + => Double -- ^ Rate parameter, also known as lambda + -> g -- ^ Generator + -> m Int +{-# INLINE poissonAtkinson #-} +poissonAtkinson lambda gen = loop + where loop :: m Int + loop = do + bigU <- uniformDouble01M gen + let x = (alpha - log ((1.0 - bigU) / bigU)) / bbeta + if x < (-0.5) + then loop + else do + bigV <- uniformDouble01M gen + let n = floor (x + 0.5) :: Int + y = alpha - bbeta * x + logFacN = logFactorial n + lhs = y + log (bigV / (1.0 + exp y)**2) + rhs = bigK + fromIntegral n * logLambda - logFacN + if lhs <= rhs + then return n + else loop + bigC,alpha,bbeta,bigK,logLambda :: Double + bigC = 0.767 - 3.36 / lambda + bbeta = pi / sqrt (3.0 * lambda) + alpha = bbeta * lambda + bigK = log bigC - lambda - log bbeta + logLambda = log lambda + -- $references -- -- * Doornik, J.A. (2005) An improved ziggurat method to generate @@ -500,3 +578,12 @@ -- 1988) 216. -- Here's an example of how the algorithm's sampling regions look -- ![Something](docs/RecreateFigure.svg) +-- +-- * Devroye, L. (1986) Non-uniform Random Variate Generation. +-- Springer Verlag. Chapter 10: Discrete Univariate Distributions. +-- +-- +-- * Robert, C.P. & Casella, G. Monte Carlo Statistical Methods. +-- Springer Texts in Statistics. Algorithm A6: Atkinson's Method +-- for Generating Poisson Random Variables. +-- diff --git a/bench/Benchmark.hs b/bench/Benchmark.hs --- a/bench/Benchmark.hs +++ b/bench/Benchmark.hs @@ -108,6 +108,10 @@ ] ] ] + , bgroup "poisson" + [ bench (show lam) $ whnfIO $ loop iter (poisson lam mwc) + | lam <- [0.1, 1, 8, 12, 100, 1000, 10000] + ] -- Test sampling performance. Table creation must be floated out! , bgroup "CT/gen" $ concat [ [ bench ("uniform "++show i) $ whnfIO $ loop iter (genFromTable tbl mwc) diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -1,9 +1,13 @@ -## Changes in 0.15.2.0 +## Changes in 0.15.3.0 [2025-12-30] + * `poisson` samples for Poisson distribution is added + +## Changes in 0.15.2.0 [2025-01-18] + * Support for `random-1.3`. -## Changes in 0.15.1.0 +## Changes in 0.15.1.0 [2024-07-09] * Additon of binomial sampler using the rejection sampling method in Kachitvichyanukul, V. and Schmeiser, B. W. Binomial Random @@ -12,17 +16,17 @@ efficient basis for e.g. the beta binomial distribution: `beta a b g >>= \p -> binomial n p g`. -## Changes in 0.15.0.2 +## Changes in 0.15.0.2 [2021-08-15] * Doctests on 32-bit platforms are fixed. (#79) -## Changes in 0.15.0.1 +## Changes in 0.15.0.1 [2020-08-08] * Bug in generation of Int/Word in both uniform and uniformR is fixed. (#75) -## Changes in 0.15.0.0 +## Changes in 0.15.0.0 [2020-07-31] * `withSystemRandomST` and `createSystemSeed` are added. diff --git a/mwc-random.cabal b/mwc-random.cabal --- a/mwc-random.cabal +++ b/mwc-random.cabal @@ -1,7 +1,7 @@ cabal-version: 3.0 build-type: Simple name: mwc-random -version: 0.15.2.0 +version: 0.15.3.0 license: BSD-2-Clause license-file: LICENSE copyright: 2009, 2010, 2011 Bryan O'Sullivan @@ -45,15 +45,15 @@ || ==9.0.2 || ==9.2.8 || ==9.4.8 - || ==9.6.6 + || ==9.6.7 || ==9.8.4 - || ==9.10.1 - || ==9.12.1 + || ==9.10.2 + || ==9.12.2 source-repository head type: git - location: git://github.com/haskell/mwc-random + location: https://github.com/haskell/mwc-random.git flag BenchPAPI Description: Enable building of benchmarks which use instruction counters. diff --git a/tests/props.hs b/tests/props.hs --- a/tests/props.hs +++ b/tests/props.hs @@ -5,7 +5,7 @@ import Data.Proxy import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Unboxed.Mutable as MVU -import Numeric.SpecFunctions (logChoose,incompleteGamma,log1p) +import Numeric.SpecFunctions (logChoose,incompleteGamma,log1p,logFactorial) import Test.Tasty import Test.Tasty.QuickCheck @@ -80,6 +80,14 @@ assertEqual "[Word32]" xs golden , testCase "beta binomial mean" $ prop_betaBinomialMean , testProperty "binomial is binomial" $ prop_binomial_PMF n_per_bin p_val g0 + , testProperty "Poisson is Poisson" $ prop_Poisson_PMF n_per_bin p_val g0 + , testCase "poisson mean and variance of 1" $ prop_poissonMeanAndVar 1 + , testCase "poisson mean and variance of 5" $ prop_poissonMeanAndVar 5 + , testCase "poisson mean and variance of 40" $ prop_poissonMeanAndVar 40 + , testCase "poisson mean and variance of 150" $ prop_poissonMeanAndVar 150 + , testCase "poisson 0" $ do + n <- poisson 0 g0 + 0 @=? n ] updateGenState :: GenIO -> IO () @@ -206,7 +214,15 @@ let x1 = fromIntegral nTrials * alpha / (alpha + delta) assertBool ("Mean is " ++ show x1 ++ " but estimated as " ++ show m) (abs (m - x1) < 0.001) - +prop_poissonMeanAndVar :: Double -> IO () +prop_poissonMeanAndVar lambda = do + gen <- create + ss <- replicateM nSamples $ poisson lambda gen + let m = fromIntegral (sum ss) / fromIntegral nSamples :: Double + let v = (fromIntegral (sum (map (^ 2) ss)) / fromIntegral nSamples) - (m ** 2) + assertBool + ("True mean and var: " ++ show lambda ++ " but estimated as " ++ show (m, v) ++ " respectively") + ((abs (lambda - m) / lambda < 0.025) && (abs (lambda - v) / lambda < 0.025)) -- Test that `binomial` really samples from binomial distribution. @@ -246,6 +262,34 @@ $ counterexample ("chi2 = " ++ show logL) $ significance > p_val +-- Similar test for Poisson distribution. We build histogram by moving +-- all values >λ+2σ into overflow bin. +prop_Poisson_PMF :: NPerBin -> PValue -> GenIO -> Property +prop_Poisson_PMF (NPerBin n_per_bin) (PValue p_val) g = property $ do + lam <- choose (0, 100.0) -- Poisson rate + let max_n = ceiling $ lam + 2 * sqrt lam -- Maximum histogram bin + n_samples = max_n * n_per_bin -- Number of samples to generate + n_samples' = fromIntegral n_samples + pure $ ioProperty $ do + hist <- do + buf <- MVU.new (max_n + 1) + replicateM_ n_samples $ do + k <- poisson lam g + MVU.modify buf (+(1::Int)) (min max_n k) + U.unsafeFreeze buf + let likelihood _ 0 + = 0 + likelihood k (fromIntegral -> n_obs) + = n_obs * (log (n_obs / n_samples') - p) + where + p | k == max_n = logComlCdfPoisson lam k + | otherwise = logProbPoisson lam k + let logL = 2 * U.sum (U.imap likelihood hist) + let significance = 1 - cumulativeChi2 max_n logL + pure $ counterexample ("lambda = " ++ show lam) + $ counterexample ("p-val = " ++ show significance) + $ counterexample ("chi2 = " ++ show logL) + $ significance > p_val ---------------------------------------------------------------- -- Statistical helpers @@ -258,6 +302,12 @@ where k' = fromIntegral k nk' = fromIntegral $ n - k + +logProbPoisson :: Double -> Int -> Double +logProbPoisson lam i = log lam * fromIntegral i - logFactorial i - lam + +logComlCdfPoisson :: Double -> Int -> Double +logComlCdfPoisson lam i = log $ incompleteGamma (fromIntegral (i :: Int)) lam cumulativeChi2 :: Int -> Double -> Double