Commits

Aleksey Khudyakov  committed 11ae934

Add table generation for poisson distribution

  • Participants
  • Parent commits 5581c03

Comments (0)

Files changed (5)

File System/Random/MWC/CondensedTable.hs

   , tableFromProbabilities
   , tableFromWeights
   , tableFromIntWeights
+    -- ** Disrete distributions
+  , tablePoisson
   ) where
 
 import Control.Arrow           (second,(***))
 import Data.Int
 import Data.Bits
 import qualified Data.Vector.Generic         as G
+import           Data.Vector.Generic           ((++),(!))
 import qualified Data.Vector.Generic.Mutable as M
 import qualified Data.Vector.Unboxed         as U
 import qualified Data.Vector                 as V
 import Data.Vector.Generic (Vector)
 
+import Prelude hiding ((++))
+
 import System.Random.MWC
 
 
   return arr
 
 
+-- | Create lookup table for poisson distibution.
+tablePoisson :: Double -> CondensedTableU Int
+tablePoisson = tableFromProbabilities . make
+  where
+    make lam
+      | lam < 0    = error "System.Random.MWC.CondesedTable.tablePoisson: negative lambda"
+      | lam < 22.8 = U.unfoldr unfoldForward (exp (-lam), 0)
+      | otherwise  = U.unfoldr unfoldForward (pMax, nMax)
+                  ++ U.tail (U.unfoldr unfoldBackward (pMax, nMax))
+      where
+        -- Number with highest probability and its probability
+        nMax = floor lam :: Int
+        pMax = let c = lam * exp( -lam / fromIntegral nMax )
+               in  U.foldl' (\p i -> p * c / i) 1 (U.enumFromN 1 nMax)
+        -- Build probability list
+        unfoldForward (p,i)
+          | p < minP  = Nothing
+          | otherwise = Just ( (i,p)
+                             , (p * lam / fromIntegral (i+1), i+1)
+                             )
+        -- Go down
+        unfoldBackward (p,i)
+          | p < minP  = Nothing
+          | otherwise = Just ( (i,p)
+                             , (p / lam * fromIntegral i, i-1)
+                             )
+    minP = 1.1641532182693481e-10 -- 2**(-33)
+
+
 -- $references
 --
 -- * Wang, J.; Tsang, W. W.; G. Marsaglia (2004), Fast Generation of
 --   Discrete Random Variables, /Journal of Statistical Software,
 --   American Statistical Association/, vol. 11(i03).
 --   <http://ideas.repec.org/a/jss/jstsof/11i03.html>
-

File benchmarks/Benchmark.hs

         , bench "gamma,a>1"   (gamma 2   1   mwc :: IO Double)
         , bench "chiSquare"   (chiSquare 4   mwc :: IO Double)
         ]
-      , bgroup "CT/gen"  $
-        map (\i -> bench ("uniform "++show i) (genFromTable (makeTableUniform i) mwc :: IO Int)) [2..10]
-      , bgroup "CT/table" $
-        map (\i -> bench ("makeTableUniform " ++ show i) (whnf makeTableUniform i)) [2..30]
+      , bgroup "CT/gen" $ concat
+        [ [ bench ("uniform "++show i)     (genFromTable (makeTableUniform i) mwc :: IO Int)
+          | i <- [2..10]
+          ]
+        , [ bench ("poisson " ++ show l)   (genFromTable (tablePoisson l) mwc :: IO Int)
+          | l <- [0.01, 0.2, 0.8, 1.3, 2.4, 8, 12, 100, 1000]
+          ]
+        ]
+      , bgroup "CT/table" $ concat
+        [ [ bench ("uniform " ++ show i) $ whnf makeTableUniform i
+          | i <- [2..30]
+          ]
+        , [ bench ("poisson " ++ show l) $ whnf tablePoisson l
+          | l <- [0.01, 0.2, 0.8, 1.3, 2.4, 8, 12, 100, 1000]
+          ]
+        ]
       ]
     , bgroup "random"
       [

File test/ChiSquare.hs

 
 import Data.Typeable
 import Data.Word
+import Data.List (find)
 import qualified Data.Vector.Unboxed              as U
 import qualified Data.Vector.Unboxed.Mutable      as M
 import qualified System.Random.MWC                as MWC
 import qualified System.Random.MWC.CondensedTable as MWC
 
 import Statistics.Test.ChiSquared
+import Statistics.Distribution
+import Statistics.Distribution.Poisson
 
 import Test.HUnit hiding (Test)
 import Test.Framework
 import Test.Framework.Providers.HUnit
 
+
+
 ----------------------------------------------------------------
 
 tests :: MWC.GenIO -> Test
   , ctableTest   [1/3 , 1/3, 1/3] g
   , ctableTest   [0.1,  0.9] g
   , ctableTest   (replicate 10 0.1) g
+    -- ** Poisson
+  , poissonTest 0.2  g
+  , poissonTest 1.32 g
+  , poissonTest 6.8  g
+  , poissonTest 100  g
   ]
 
 ----------------------------------------------------------------
               , probabilites = U.fromList ps
               }
     r <- sampleTest gen (10^4) g
+    assertEqual "Significant!" NotSignificant r
+
+-- | Test for condensed table for poissson distribution
+poissonTest :: Double -> MWC.GenIO -> Test
+poissonTest lam g =
+  testCase ("poissonTest: " ++ show lam) $ do
+    let pois      = poisson lam
+        Just nMax = find (\n -> probability pois n < 2**(-33)) [floor lam ..]
+    let gen = Generator
+              { generator    = MWC.genFromTable (MWC.tablePoisson lam)
+              , probabilites = U.generate nMax (probability pois)
+              }
+    r <- sampleTest gen (10^4) g
     assertEqual "Significant!" NotSignificant r

File test/visual.R

 
 
 view.dumps <- function() {
+  # Load random data from dist
   load.d <- function(name) read.table(name)[,1]
+  # Plots for continous distribution
   plot.d <- function(name, dens, rng) {
     smp <- load.d( name )
     plot( density(smp), xlim=rng, main=name, col='blue', lwd=2)
     hist( smp, probability=TRUE, breaks=100, add=TRUE)
     plot( dens, xlim=rng, col='red', add=TRUE, lwd=2)
   }
+  # plots for discrete distribution
+  plot.ds <- function( name, xs, prob) {
+    smp <- load.d( name )
+    h   <- hist( smp, breaks = xs - 0.5, freq=FALSE )
+    dh  <- sqrt( h$count ) / max( 1, sum( h$count ) )
+    arrows( xs, h$density + dh,
+            xs, h$density - dh,
+            angle=90, code=3, length=0.2 )
+    points( xs, prob(xs), pch='0', col='red', type='b')
+  }
   ################################################################
   # Normal
   plot.d ("distr/normal-0-1",
           function(x) dexp(x,3),
           c(-0.5, 3) )
   readline()
+  ################################################################
+  # Poisson
+  plot.ds( "distr/poisson-0.1", 0:6, function(x) dpois(x, lambda=0.1) )
+  readline()
+  #
+  plot.ds( "distr/poisson-1.0", 0:10, function(x) dpois(x, lambda=1.0) )
+  readline()
+  #
+  plot.ds( "distr/poisson-4.5", 0:20, function(x) dpois(x, lambda=4.5) )
+  readline()
+  #
+  plot.ds( "distr/poisson-30", 0:100, function(x) dpois(x, lambda=30) )
+  readline()
+  #
 }

File test/visual.hs

 import System.Directory  (createDirectoryIfMissing,setCurrentDirectory)
 import System.IO
 
-import qualified System.Random.MWC               as MWC
-import qualified System.Random.MWC.Distributions as MWC
+import qualified System.Random.MWC                as MWC
+import qualified System.Random.MWC.Distributions  as MWC
+import qualified System.Random.MWC.CondensedTable as MWC
 
 
 dumpSample :: Show a => Int -> FilePath -> IO a -> IO ()
   -- Exponential
   dumpSample n "exponential-1" $ MWC.exponential 1 g
   dumpSample n "exponential-3" $ MWC.exponential 3 g
+  -- Poisson
+  dumpSample n "poisson-0.1"   $ MWC.genFromTable (MWC.tablePoisson 0.1) g
+  dumpSample n "poisson-1.0"   $ MWC.genFromTable (MWC.tablePoisson 1.0) g
+  dumpSample n "poisson-4.5"   $ MWC.genFromTable (MWC.tablePoisson 4.5) g
+  dumpSample n "poisson-30"    $ MWC.genFromTable (MWC.tablePoisson 30)  g