Commits

Aleksey Khudyakov committed f147b95

Add table generation forn binomial distribution

  • Participants
  • Parent commits 11ae934

Comments (0)

Files changed (5)

File System/Random/MWC/CondensedTable.hs

   , tableFromIntWeights
     -- ** Disrete distributions
   , tablePoisson
+  , tableBinomial
   ) where
 
 import Control.Arrow           (second,(***))
                              )
     minP = 1.1641532182693481e-10 -- 2**(-33)
 
+-- | Create lookup table for binomial distribution
+tableBinomial :: Int            -- ^ Number of tries
+              -> Double         -- ^ Probability of success
+              -> CondensedTableU Int
+tableBinomial n p = tableFromProbabilities makeBinom
+  where 
+  makeBinom
+    | n <= 0         = error "System.Random.MWC.CondesedTable.tableBinomial: nonpositive number of tryes"
+    | p == 0         = U.singleton (0,1)
+    | p == 1         = U.singleton (n,1)
+    | p > 0 && p < 1 = U.unfoldrN (n + 1) unfolder ((1-p)^n, 0)
+    | otherwise      = error "System.Random.MWC.CondesedTable.tableBinomial: probability is out of range"
+    where
+      h = p / (1 - p)
+      unfolder (t,i) = Just ( (i,t)
+                            , (t * (fromIntegral $ n + 1 - i1) * h / fromIntegral i1, i1) )
+        where i1 = i + 1
+
 
 -- $references
 --

File benchmarks/Benchmark.hs

         , [ 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]
           ]
+        , [ bench ("binomial " ++ show p ++ " " ++ show n) (genFromTable (tableBinomial n p) mwc :: IO Int)
+          | (n,p) <- [ (4, 0.5), (10,0.1), (10,0.6), (10, 0.8), (100,0.4)]
+          ]
         ]
       , bgroup "CT/table" $ concat
         [ [ bench ("uniform " ++ show i) $ whnf makeTableUniform i
         , [ bench ("poisson " ++ show l) $ whnf tablePoisson l
           | l <- [0.01, 0.2, 0.8, 1.3, 2.4, 8, 12, 100, 1000]
           ]
+        , [ bench ("binomial " ++ show p ++ " " ++ show n) $ whnf (tableBinomial n) p
+          | (n,p) <- [ (4, 0.5), (10,0.1), (10,0.6), (10, 0.8), (100,0.4)]
+          ]
         ]
       ]
     , bgroup "random"

File test/ChiSquare.hs

 import Statistics.Test.ChiSquared
 import Statistics.Distribution
 import Statistics.Distribution.Poisson
+import Statistics.Distribution.Binomial
 
 import Test.HUnit hiding (Test)
 import Test.Framework
   , poissonTest 1.32 g
   , poissonTest 6.8  g
   , poissonTest 100  g
+    -- ** Binomial
+  , binomialTest 4   0.5 g
+  , binomialTest 10  0.1 g
+  , binomialTest 10  0.6 g
+  , binomialTest 10  0.8 g
+  , binomialTest 100 0.3 g
   ]
 
 ----------------------------------------------------------------
               , probabilites = U.generate nMax (probability pois)
               }
     r <- sampleTest gen (10^4) g
-    assertEqual "Significant!" NotSignificant r
+    assertEqual "Significant!" NotSignificant r
+    
+binomialTest :: Int -> Double -> MWC.GenIO -> Test
+binomialTest n p g =
+  testCase ("binomialTest: " ++ show p ++ " " ++ show n) $ do
+    let binom = binomial n p
+        gen = Generator
+              { generator    = MWC.genFromTable (MWC.tableBinomial n p)
+              , probabilites = U.generate (n+1) (probability binom)
+              }
+    r <- sampleTest gen (10^4) g
+    assertEqual "Significant!" NotSignificant r
+

File test/visual.R

   # plots for discrete distribution
   plot.ds <- function( name, xs, prob) {
     smp <- load.d( name )
-    h   <- hist( smp, breaks = xs - 0.5, freq=FALSE )
+    h   <- hist( smp,
+                 breaks = c( max(xs) + 0.5, xs - 0.5),
+                 freq=FALSE, main = name
+                )
     dh  <- sqrt( h$count ) / max( 1, sum( h$count ) )
     arrows( xs, h$density + dh,
             xs, h$density - dh,
   plot.ds( "distr/poisson-30", 0:100, function(x) dpois(x, lambda=30) )
   readline()
   #
+  ################################################################
+  # Binomial
+  plot.ds( "distr/binom-4-0.5", 0:4, function(x) dbinom(x, 4, 0.5) )
+  readline()
+  #
+  plot.ds( "distr/binom-10-0.1", 0:10, function(x) dbinom(x, 10, 0.1) )
+  readline()
+  #
+  plot.ds( "distr/binom-10-0.6", 0:10, function(x) dbinom(x, 10, 0.6) )
+  readline()
+  #
+  plot.ds( "distr/binom-10-0.8", 0:10, function(x) dbinom(x, 10, 0.8) )
+  readline()
+  #
 }

File test/visual.hs

 dumpSample n fname gen =
   withFile fname WriteMode $ \h -> 
     replicateM_ n (hPutStrLn h . show =<< gen)
-  
+
 main :: IO ()
 main = MWC.withSystemRandom $ \g -> do
-  let n   = 10000
+  let n   = 30000
       dir = "distr"
   createDirectoryIfMissing True dir
   setCurrentDirectory           dir
   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
+  -- Binomial
+  dumpSample n "binom-4-0.5"   $ MWC.genFromTable (MWC.tableBinomial 4  0.5) g
+  dumpSample n "binom-10-0.1"  $ MWC.genFromTable (MWC.tableBinomial 10 0.1) g  
+  dumpSample n "binom-10-0.6"  $ MWC.genFromTable (MWC.tableBinomial 10 0.6) g
+  dumpSample n "binom-10-0.8"  $ MWC.genFromTable (MWC.tableBinomial 10 0.8) g