random-hypergeometric (empty) → 0.1.0.0
raw patch · 6 files changed
+223/−0 lines, 6 filesdep +Cabaldep +QuickCheckdep +basesetup-changed
Dependencies added: Cabal, QuickCheck, base, cabal-test-quickcheck, math-functions, mwc-random, random-fu, vector
Files
- LICENSE +22/−0
- Setup.hs +2/−0
- random-hypergeometric.cabal +41/−0
- src/Data/Random/Distribution/Hypergeometric.hs +56/−0
- src/Data/Random/Distribution/Hypergeometric/Impl.hs +43/−0
- src/Data/Random/Distribution/Hypergeometric/Test.hs +59/−0
+ LICENSE view
@@ -0,0 +1,22 @@+The MIT License (MIT)++Copyright (c) 2015 Sam Rijs++Permission is hereby granted, free of charge, to any person obtaining a copy+of this software and associated documentation files (the "Software"), to deal+in the Software without restriction, including without limitation the rights+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell+copies of the Software, and to permit persons to whom the Software is+furnished to do so, subject to the following conditions:++The above copyright notice and this permission notice shall be included in all+copies or substantial portions of the Software.++THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE+SOFTWARE.+
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ random-hypergeometric.cabal view
@@ -0,0 +1,41 @@+name: random-hypergeometric+version: 0.1.0.0+synopsis: Random variate generation from hypergeometric distributions+description: The Hypergeometric distribution. This is the discrete probability+ distribution that measures the probability of /k/ successes in /l/+ trials, without replacement, from a finite population.+homepage: https://github.com/srijs/random-hypergeometric+license: MIT+license-file: LICENSE+author: Sam Rijs+maintainer: srijs@airpost.net+copyright: 2015 Sam Rijs+ 2005 Robert Kern+ 1998 Ivan Frohne+category: Math+build-type: Simple+cabal-version: >=1.10++library+ exposed-modules: Data.Random.Distribution.Hypergeometric+ other-modules: Data.Random.Distribution.Hypergeometric.Impl+ build-depends: base >=4.7 && <4.8,+ random-fu >=0.2 && <0.3,+ math-functions >=0.1 && <0.2+ hs-source-dirs: src+ default-language: Haskell2010++test-suite test+ type: detailed-0.9+ test-module: Data.Random.Distribution.Hypergeometric.Test+ other-modules: Data.Random.Distribution.Hypergeometric,+ Data.Random.Distribution.Hypergeometric.Impl+ build-depends: base >=4.7 && <4.8,+ Cabal >=1.10,+ random-fu >=0.2 && <0.3,+ math-functions >=0.1 && <0.2,+ mwc-random >=0.13 && <0.14,+ vector >=0.10 && <0.11,+ QuickCheck >=2.7 && <2.8,+ cabal-test-quickcheck >=0.1 && <0.2+ hs-source-dirs: src
+ src/Data/Random/Distribution/Hypergeometric.hs view
@@ -0,0 +1,56 @@+{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE+ MultiParamTypeClasses,+ FlexibleContexts,+ FlexibleInstances+ #-}+-- |+-- Module : Statistics.Distribution.Hypergeometric.GenVar+-- Copyright : (c) 2015 Sam Rijs,+-- (c) 2005 Robert Kern,+-- (c) 1998 Ivan Frohne+-- License : MIT+--+-- Maintainer : srijs@airpost.net+-- Stability : experimental+--+-- The parameters of the distribution describe /k/ elements chosen+-- from a population of /l/, with /m/ elements of one type, and+-- /l/-/m/ of the other (all are positive integers).++module Data.Random.Distribution.Hypergeometric+ ( Hypergeometric+ -- ** Constructors+ , hypergeometric+ -- ** Accessors+ , getM, getL, getK+ -- ** Variate Generation+ , hypergeometricVar+ , hypergeometricVarT+ ) where++import Data.Random.RVar+import Data.Random.Distribution+import Data.Random.Distribution.Uniform++import Data.Random.Distribution.Hypergeometric.Impl++data Hypergeometric t = Hypergeometric { getK :: !t, getL :: !t, getM :: !t }++-- | Constructs a hypergeometric distribution from the parameters /k/, /l/ and /m/.+-- Fails if /l/ is negative, /k/ is not in [0,/l/] or /m/ is not in [0,/l/].+hypergeometric :: (Num a, Ord a) => a -> a -> a -> Hypergeometric a+hypergeometric k l m + | l < 0 = error "l must not be negative"+ | m < 0 || m > l = error "m must be in [0,l]"+ | k < 0 || k > l = error "k must be in [0,l]"+ | otherwise = Hypergeometric k l m++hypergeometricVar :: (Num a, Ord a, Distribution Hypergeometric a) => a -> a -> a -> RVar a+hypergeometricVar = hypergeometricVarT++hypergeometricVarT :: (Num a, Ord a, Distribution Hypergeometric a) => a -> a -> a -> RVarT m a+hypergeometricVarT k l m = rvarT (hypergeometric k l m)++instance (Integral t) => Distribution Hypergeometric t where+ rvarT (Hypergeometric k l m) = rhyper (k, l, m)
+ src/Data/Random/Distribution/Hypergeometric/Impl.hs view
@@ -0,0 +1,43 @@+module Data.Random.Distribution.Hypergeometric.Impl where++import Data.Random.RVar+import Data.Random.Distribution+import Data.Random.Distribution.Uniform++import Numeric.SpecFunctions (logGamma)++d1 = 1.7155277699214135 -- 2*sqrt(2/e)+d2 = 0.8989161620588988 -- 3 - 2*sqrt(3/e)++rhyper :: (Integral a) => (a, a, a) -> RVarT m a+rhyper (sample, popsize, good) = return . fix2 . fix1 =<< loop+ where fix1 z = if good > bad then m - z else z+ fix2 z = if m < sample then good - z else z+ mingoodbad = min good bad+ bad = popsize - good+ maxgoodbad = max good bad+ m = min sample (popsize - sample)+ gamma z = sum $ map (logGamma . fromIntegral)+ [ z + 1, mingoodbad - z + 1+ , m - z + 1, maxgoodbad - m + z + 1 ]+ d4 = fromIntegral mingoodbad / fromIntegral popsize+ d5 = 1 - d4+ d6 = fromIntegral m * d4 + 0.5+ d7 = sqrt $ fromIntegral (popsize - m) * fromIntegral sample * d4 * d5 / fromIntegral (popsize - 1) + 0.5+ d8 = d1 * d7 + d2+ d9 = floor $ fromIntegral (m + 1) * fromIntegral (mingoodbad + 1) / fromIntegral (popsize + 2)+ d10 = gamma d9+ d11 = fromIntegral $ min (min m mingoodbad + 1) (floor (d6 + 16 * d7))+ -- 16 for 16-decimal-digit precision in d1 and d2+ loop = do+ x <- doubleUniform 0 1+ y <- doubleUniform 0 1+ case () of + _ | w < 0 || w >= d11 -> loop -- fast rejection+ | x * (4 - x) - 3 <= t -> return z -- fast acceptance+ | x * (x - t) >= 1 -> loop -- fast rejection+ | 2 * (log x) <= t -> return z -- acceptance+ | otherwise -> loop -- rejection+ where w = d6 + d8 * (y - 0.5) / x+ z = floor w+ t = d10 - gamma z
+ src/Data/Random/Distribution/Hypergeometric/Test.hs view
@@ -0,0 +1,59 @@+module Data.Random.Distribution.Hypergeometric.Test where++import Data.Int++import Data.Random+import Data.Random.Sample (sampleFrom)+import Data.Random.Distribution.Hypergeometric++import System.Random.MWC++import Data.Vector (fromList)++import Test.QuickCheck+import Test.QuickCheck.Monadic++import Distribution.TestSuite.QuickCheck++instance Arbitrary Seed where+ arbitrary = return . toSeed . fromList =<< vector 258++testArbitraryRandom :: (GenIO -> IO Bool) -> Seed -> Property+testArbitraryRandom f s = monadicIO $ assert =<< run (f =<< restore s)++tests :: IO [Test]+tests = return++ [ testGroup "identities"++ [ testGroup "small"++ [ testProperty "draw nil" $ \p -> testArbitraryRandom $ \g -> do+ let i = (getSmall . getPositive) p :: Int64+ v <- sampleFrom g $ hypergeometric i i 0+ return $ v == 0++ , testProperty "draw all" $ \p -> testArbitraryRandom $ \g -> do+ let i = (getSmall . getPositive) p :: Int64+ v <- sampleFrom g $ hypergeometric i i i+ return $ v == i++ ]++ , testGroup "large"++ [ testProperty "draw nil" $ \p -> testArbitraryRandom $ \g -> do+ let i = (getLarge . getPositive) p :: Int64+ v <- sampleFrom g $ hypergeometric i i 0+ return $ v == 0++ , testProperty "draw all" $ \p -> testArbitraryRandom $ \g -> do+ let i = (getLarge . getPositive) p :: Int64+ v <- sampleFrom g $ hypergeometric i i i+ return $ v == i++ ]++ ]++ ]