Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse code

Add table generation forn binomial distribution

  • Loading branch information...
commit c47610bd3b22804588bfa84720c2bb4391ac7404 1 parent e62ae0d
Aleksey Khudyakov authored February 09, 2012
19  System/Random/MWC/CondensedTable.hs
@@ -17,6 +17,7 @@ module System.Random.MWC.CondensedTable (
17 17
   , tableFromIntWeights
18 18
     -- ** Disrete distributions
19 19
   , tablePoisson
  20
+  , tableBinomial
20 21
   ) where
21 22
 
22 23
 import Control.Arrow           (second,(***))
@@ -226,6 +227,24 @@ tablePoisson = tableFromProbabilities . make
226 227
                              )
227 228
     minP = 1.1641532182693481e-10 -- 2**(-33)
228 229
 
  230
+-- | Create lookup table for binomial distribution
  231
+tableBinomial :: Int            -- ^ Number of tries
  232
+              -> Double         -- ^ Probability of success
  233
+              -> CondensedTableU Int
  234
+tableBinomial n p = tableFromProbabilities makeBinom
  235
+  where 
  236
+  makeBinom
  237
+    | n <= 0         = error "System.Random.MWC.CondesedTable.tableBinomial: nonpositive number of tryes"
  238
+    | p == 0         = U.singleton (0,1)
  239
+    | p == 1         = U.singleton (n,1)
  240
+    | p > 0 && p < 1 = U.unfoldrN (n + 1) unfolder ((1-p)^n, 0)
  241
+    | otherwise      = error "System.Random.MWC.CondesedTable.tableBinomial: probability is out of range"
  242
+    where
  243
+      h = p / (1 - p)
  244
+      unfolder (t,i) = Just ( (i,t)
  245
+                            , (t * (fromIntegral $ n + 1 - i1) * h / fromIntegral i1, i1) )
  246
+        where i1 = i + 1
  247
+
229 248
 
230 249
 -- $references
231 250
 --
6  benchmarks/Benchmark.hs
@@ -69,6 +69,9 @@ main = do
69 69
         , [ bench ("poisson " ++ show l)   (genFromTable (tablePoisson l) mwc :: IO Int)
70 70
           | l <- [0.01, 0.2, 0.8, 1.3, 2.4, 8, 12, 100, 1000]
71 71
           ]
  72
+        , [ bench ("binomial " ++ show p ++ " " ++ show n) (genFromTable (tableBinomial n p) mwc :: IO Int)
  73
+          | (n,p) <- [ (4, 0.5), (10,0.1), (10,0.6), (10, 0.8), (100,0.4)]
  74
+          ]
72 75
         ]
73 76
       , bgroup "CT/table" $ concat
74 77
         [ [ bench ("uniform " ++ show i) $ whnf makeTableUniform i
@@ -77,6 +80,9 @@ main = do
77 80
         , [ bench ("poisson " ++ show l) $ whnf tablePoisson l
78 81
           | l <- [0.01, 0.2, 0.8, 1.3, 2.4, 8, 12, 100, 1000]
79 82
           ]
  83
+        , [ bench ("binomial " ++ show p ++ " " ++ show n) $ whnf (tableBinomial n) p
  84
+          | (n,p) <- [ (4, 0.5), (10,0.1), (10,0.6), (10, 0.8), (100,0.4)]
  85
+          ]
80 86
         ]
81 87
       ]
82 88
     , bgroup "random"
21  test/ChiSquare.hs
@@ -18,6 +18,7 @@ import qualified System.Random.MWC.CondensedTable as MWC
18 18
 import Statistics.Test.ChiSquared
19 19
 import Statistics.Distribution
20 20
 import Statistics.Distribution.Poisson
  21
+import Statistics.Distribution.Binomial
21 22
 
22 23
 import Test.HUnit hiding (Test)
23 24
 import Test.Framework
@@ -48,6 +49,12 @@ tests g = testGroup "Chi squared tests"
48 49
   , poissonTest 1.32 g
49 50
   , poissonTest 6.8  g
50 51
   , poissonTest 100  g
  52
+    -- ** Binomial
  53
+  , binomialTest 4   0.5 g
  54
+  , binomialTest 10  0.1 g
  55
+  , binomialTest 10  0.6 g
  56
+  , binomialTest 10  0.8 g
  57
+  , binomialTest 100 0.3 g
51 58
   ]
52 59
 
53 60
 ----------------------------------------------------------------
@@ -117,4 +124,16 @@ poissonTest lam g =
117 124
               , probabilites = U.generate nMax (probability pois)
118 125
               }
119 126
     r <- sampleTest gen (10^4) g
120  
-    assertEqual "Significant!" NotSignificant r
  127
+    assertEqual "Significant!" NotSignificant r
  128
+    
  129
+binomialTest :: Int -> Double -> MWC.GenIO -> Test
  130
+binomialTest n p g =
  131
+  testCase ("binomialTest: " ++ show p ++ " " ++ show n) $ do
  132
+    let binom = binomial n p
  133
+        gen = Generator
  134
+              { generator    = MWC.genFromTable (MWC.tableBinomial n p)
  135
+              , probabilites = U.generate (n+1) (probability binom)
  136
+              }
  137
+    r <- sampleTest gen (10^4) g
  138
+    assertEqual "Significant!" NotSignificant r
  139
+
19  test/visual.R
@@ -15,7 +15,10 @@ view.dumps <- function() {
15 15
   # plots for discrete distribution
16 16
   plot.ds <- function( name, xs, prob) {
17 17
     smp <- load.d( name )
18  
-    h   <- hist( smp, breaks = xs - 0.5, freq=FALSE )
  18
+    h   <- hist( smp,
  19
+                 breaks = c( max(xs) + 0.5, xs - 0.5),
  20
+                 freq=FALSE, main = name
  21
+                )
19 22
     dh  <- sqrt( h$count ) / max( 1, sum( h$count ) )
20 23
     arrows( xs, h$density + dh,
21 24
             xs, h$density - dh,
@@ -85,4 +88,18 @@ view.dumps <- function() {
85 88
   plot.ds( "distr/poisson-30", 0:100, function(x) dpois(x, lambda=30) )
86 89
   readline()
87 90
   #
  91
+  ################################################################
  92
+  # Binomial
  93
+  plot.ds( "distr/binom-4-0.5", 0:4, function(x) dbinom(x, 4, 0.5) )
  94
+  readline()
  95
+  #
  96
+  plot.ds( "distr/binom-10-0.1", 0:10, function(x) dbinom(x, 10, 0.1) )
  97
+  readline()
  98
+  #
  99
+  plot.ds( "distr/binom-10-0.6", 0:10, function(x) dbinom(x, 10, 0.6) )
  100
+  readline()
  101
+  #
  102
+  plot.ds( "distr/binom-10-0.8", 0:10, function(x) dbinom(x, 10, 0.8) )
  103
+  readline()
  104
+  #
88 105
 }
9  test/visual.hs
@@ -13,10 +13,10 @@ dumpSample :: Show a => Int -> FilePath -> IO a -> IO ()
13 13
 dumpSample n fname gen =
14 14
   withFile fname WriteMode $ \h -> 
15 15
     replicateM_ n (hPutStrLn h . show =<< gen)
16  
-  
  16
+
17 17
 main :: IO ()
18 18
 main = MWC.withSystemRandom $ \g -> do
19  
-  let n   = 10000
  19
+  let n   = 30000
20 20
       dir = "distr"
21 21
   createDirectoryIfMissing True dir
22 22
   setCurrentDirectory           dir
@@ -37,3 +37,8 @@ main = MWC.withSystemRandom $ \g -> do
37 37
   dumpSample n "poisson-1.0"   $ MWC.genFromTable (MWC.tablePoisson 1.0) g
38 38
   dumpSample n "poisson-4.5"   $ MWC.genFromTable (MWC.tablePoisson 4.5) g
39 39
   dumpSample n "poisson-30"    $ MWC.genFromTable (MWC.tablePoisson 30)  g
  40
+  -- Binomial
  41
+  dumpSample n "binom-4-0.5"   $ MWC.genFromTable (MWC.tableBinomial 4  0.5) g
  42
+  dumpSample n "binom-10-0.1"  $ MWC.genFromTable (MWC.tableBinomial 10 0.1) g  
  43
+  dumpSample n "binom-10-0.6"  $ MWC.genFromTable (MWC.tableBinomial 10 0.6) g
  44
+  dumpSample n "binom-10-0.8"  $ MWC.genFromTable (MWC.tableBinomial 10 0.8) g

0 notes on commit c47610b

Please sign in to comment.
Something went wrong with that request. Please try again.