statistics-hypergeometric-genvar (empty) → 0.1.0.0

raw patch · 4 files changed

+120/−0 lines, 4 filesdep +basedep +math-functionsdep +mwc-randomsetup-changed

Dependencies added: base, math-functions, mwc-random, primitive, statistics

Files

• LICENSE +20/−0
• Setup.hs +2/−0
• src/Statistics/Distribution/Hypergeometric/GenVar.hs +68/−0
• statistics-hypergeometric-genvar.cabal +30/−0

@@ -0,0 +1,20 @@+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.

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

@@ -0,0 +1,68 @@+{-# LANGUAGE Trustworthy #-}+-- |+-- 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+--++module Statistics.Distribution.Hypergeometric.GenVar (genVar) where++import Control.Monad.Primitive (PrimMonad, PrimState)++import Statistics.Distribution (ContGen(..), DiscreteGen(..))+import Statistics.Distribution.Hypergeometric (HypergeometricDistribution, hdM, hdL, hdK)++import Numeric.SpecFunctions (logGamma)++import System.Random.MWC (Gen, uniform)++d1 = 1.7155277699214135 -- 2*sqrt(2/e)+d2 = 0.8989161620588988 -- 3 - 2*sqrt(3/e)++instance ContGen HypergeometricDistribution+  where genContVar d g = return . fromIntegral =<< genVar (hdM d, hdL d, hdK d) g++instance DiscreteGen HypergeometricDistribution+  where genDiscreteVar d g = genVar (hdM d, hdL d, hdK d) g ++-- | Random variates from the hypergeometric distribution /(m, l, k)/.+-- Returns the number of white balls drawn when /k/ balls+-- are drawn at random from an urn containing /m/ white+-- and /l-m/ black balls.+genVar :: (PrimMonad m, Integral a) => (a, a, a) -> Gen (PrimState m) -> m a+genVar (good, popsize, sample) gen = 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 <- uniform gen+          y <- uniform gen+          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

@@ -0,0 +1,30 @@+name:                statistics-hypergeometric-genvar+version:             0.1.0.0+synopsis:            Random variate generation from hypergeometric distributions+description:         Supplements the @Statistics.Distribution.Hypergeometric@+                     module from the excellent @statistics@ package.+                     .+                     Implements the interfaces @ContGen@ and @DiscreteGen@ for+                     @HypergeometricDistribution@. Provides it's own function+                     @genVar@ generalised over @Integral@.+license:             MIT+license-file:        LICENSE+author:              Sam Rijs <srijs@airpost.net>+maintainer:          Sam Rijs <srijs@airpost.net>+homepage:            https://github.com/srijs/statistics-hypergeometric-genvar+copyright:           2015 Sam Rijs+                     2005 Robert Kern+                     1998 Ivan Frohne+category:            Math, Statistics+build-type:          Simple+cabal-version:       >=1.10++library+  exposed-modules:     Statistics.Distribution.Hypergeometric.GenVar+  build-depends:       base           >=4.7 && <4.8,+                       primitive      >=0.5 && <0.6,+                       mwc-random     >=0.13 && <0.14,+                       math-functions >=0.1 && <0.2,+                       statistics     >=0.13 && <0.14+  hs-source-dirs:      src+  default-language:    Haskell2010