Aleksey Khudyakov avatar Aleksey Khudyakov committed 4b876ac

Add tests to test generators
* QC test to check whether all values produced by uniformR lie
in range
* Chi squared test to test uniformity of distribution of uniformR
It's not perfect but should work if underlieing generator has
uniform distribution

Comments (0)

Files changed (2)

+-- QC tests for random number generators
+--
+-- Require QuickCheck >= 2.2
+
+import Control.Applicative
+import Test.QuickCheck
+import Test.QuickCheck.Arbitrary
+import Test.QuickCheck.Monadic
+import System.Random.MWC
+
+import Data.Word (Word8,Word16,Word32,Word64,Word)
+import Data.Int  (Int8, Int16, Int32, Int64 )
+
+
+
+-- Ordered pair (x,y) for which x <= y
+newtype OrderedPair a = OrderedPair (a,a)
+                        deriving Show
+instance (Ord a, Arbitrary a) => Arbitrary (OrderedPair a) where
+  arbitrary = OrderedPair <$> suchThat arbitrary (\(x,y) -> x <= y)
+
+
+----------------------------------------------------------------
+-- Test that values generated with uniformR never lie outside range.
+
+prop_InRange :: (Variate a, Ord a) => GenIO -> OrderedPair a -> Property
+prop_InRange g (OrderedPair (x1,x2)) = monadicIO $ do
+  r <- run $ uniformR (x1,x2) g
+  assert (x1 <= r && r <= x2)
+
+type InRange a = OrderedPair a -> Property
+
+test_InRange :: IO ()
+test_InRange = do
+  g <- create
+  let q :: (Testable prop) => prop -> IO ()
+      q = quickCheck
+  putStrLn "Int8"   >> q (prop_InRange g :: InRange Int8)
+  putStrLn "Int16"  >> q (prop_InRange g :: InRange Int16)
+  putStrLn "Int32"  >> q (prop_InRange g :: InRange Int32)
+  putStrLn "Int64"  >> q (prop_InRange g :: InRange Int64)
+  putStrLn "Word8"  >> q (prop_InRange g :: InRange Word8)
+  putStrLn "Word16" >> q (prop_InRange g :: InRange Word16)
+  putStrLn "Word32" >> q (prop_InRange g :: InRange Word32)
+  putStrLn "Word64" >> q (prop_InRange g :: InRange Word64)
+  putStrLn "Int"    >> q (prop_InRange g :: InRange Int)
+  putStrLn "Word64" >> q (prop_InRange g :: InRange Word)
+  putStrLn "Float"  >> q (prop_InRange g :: InRange Float)
+  putStrLn "Double" >> q (prop_InRange g :: InRange Double)
+-- Tests for testing uniformity of distributions
+--
+-- Require statistics >= 0.7
+
+{-# LANGUAGE BangPatterns #-}
+
+import Data.Word
+import Data.Int
+import qualified Data.Vector.Generic         as G
+import qualified Data.Vector.Unboxed         as U
+import qualified Data.Vector.Unboxed.Mutable as M
+
+import Statistics.Distribution
+import Statistics.Distribution.ChiSquared
+
+import System.Random.MWC
+
+import Text.Printf
+
+
+-- Fill vector with number of occurences of random number in range.
+-- Uses uniformR for generation of random numbers
+fill :: (Variate a, U.Unbox a, Integral a) => 
+        Int                     -- Expected number of items in each bin
+     -> (a,a)                   -- Range for values
+     -> GenIO                   -- Generator
+     -> IO (U.Vector Int)
+fill n rng@(x1,x2) g = do
+  let l = fromIntegral x2 - fromIntegral x1 + 1
+  v <- M.newWith l 0
+  let loop k | k == n*l  = return ()
+             | otherwise = do x <- fromIntegral `fmap` uniformR rng g
+                              M.write v x . (+1) =<< M.read v x
+                              loop (k+1)
+  loop 0
+  G.unsafeFreeze v
+
+-- Calculate χ² statistics for vector of number occurences for
+-- hypotheshys that each bin has equal probability
+chi2uniform :: U.Vector Int -> Double
+chi2uniform v = (U.sum $ U.map (sqr . subtract μ . fromIntegral) v) / μ 
+  where
+    n   = U.length v
+    tot = U.sum v
+    μ   = fromIntegral tot / fromIntegral n
+    sqr x = x * x
+
+-- Perform χ² on vector of number of occurences
+checkChi2 :: Double             -- Desired significance level
+          -> U.Vector Int       -- Vector of values
+          -> IO ()
+checkChi2 p v = do
+  let χ2   = chi2uniform v      -- Observed χ²
+      ndf  = U.length v - 1     -- N degrees of freedom
+      d    = chiSquared ndf     -- Theoretical distribution
+      pLow = cumulative d χ2
+      pHi  = 1 - pLow
+  
+  putStrLn $ if pLow > p && (1-pLow) > p then "OK" else "* FAILED *"
+  printf "  significance = %.3f\n"   p
+  printf "  χ²/ndf = %.3f\n"        (χ2 / fromIntegral ndf)
+  printf "  p(χ² < observed) = %.3g\n" pLow
+  printf "  p(χ² > observed) = %.3g\n" pHi
+
+
+main :: IO ()
+main = do
+  putStrLn "(0,255) Word8"
+  v1 <- withSystemRandom $ fill 10000 (0,255 :: Word8)
+  checkChi2 0.05 v1
+  checkChi2 0.01 v1
+  putStrLn ""
+  ----------------------------------------
+  putStrLn "(0,254) Word8"
+  v1 <- withSystemRandom $ fill 10000 (0,254 :: Word8)
+  checkChi2 0.05 v1
+  checkChi2 0.01 v1
+  putStrLn ""
+  ----------------------------------------
+  putStrLn "(0,10) Word8"
+  v1 <- withSystemRandom $ fill 10000 (0,10 :: Word8)
+  checkChi2 0.05 v1
+  checkChi2 0.01 v1
+  putStrLn ""
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.