mwc-random 0.15.2.0 → 0.15.3.0

raw patch · 5 files changed

+158/−13 lines, 5 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• System/Random/MWC/Distributions.hs +88/−1
• bench/Benchmark.hs +4/−0
• changelog.md +9/−5
• mwc-random.cabal +5/−5
• tests/props.hs +52/−2

@@ -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. <https://dl.acm.org/doi/pdf/10.1145/42372.42381> --   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.+--   <http://https://luc.devroye.org/chapter_ten.pdf>+-- +-- * Robert, C.P. & Casella, G. Monte Carlo Statistical Methods.+--   Springer Texts in Statistics. Algorithm A6: Atkinson's Method+--   for Generating Poisson Random Variables.+--   <https://mcube.lab.nycu.edu.tw/~cfung/docs/books/robert2004monte_carlo_statistical_methods.pdf>

@@ -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)

@@ -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. 

@@ -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.

@@ -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