approx-rand-test 0.0.4 → 0.1.0

raw patch · 9 files changed

+467/−129 lines, 9 filesdep +Chartdep +colourdep +containersdep −resourcetdep ~conduitdep ~statisticsdep ~vectorPVP ok

version bump matches the API change (PVP)

Dependencies added: Chart, colour, containers, data-accessor, filepath, vector-algorithms

Dependencies removed: resourcet

Dependency ranges changed: conduit, statistics, vector

API changes (from Hackage documentation)

Files

• approx-rand-test.cabal +35/−15
• src/Statistics/Test/ApproxRand.hs +69/−55
• src/Statistics/Test/Correction.hs +42/−0
• utils/Histogram.hs +33/−0
• utils/TextHistogram.hs +30/−0
• utils/approx-rand-test-paired.hs +74/−31
• utils/approx-rand-test.hs +65/−28
• utils/cairo/CairoHistogram.hs +104/−0
• utils/nocairo/CairoHistogram.hs +15/−0

@@ -1,5 +1,5 @@ Name:          approx-rand-test-Version:       0.0.4+Version:       0.1.0 License:       OtherLicense License-file:  LICENSE Copyright:     Copyright 2012 Daniël de Kok@@ -11,6 +11,9 @@ Description:   Utility to perform approximate randomization tests. Cabal-Version: >= 1.8 Build-Type:    Simple+Extra-Source-Files:+  utils/cairo/CairoHistogram.hs+  utils/nocairo/CairoHistogram.hs  Source-Repository head   Type:     git@@ -19,45 +22,62 @@ Source-Repository this   Type:     git   Location: git://github.com/danieldk/approx-rand-test.git-  Tag:      0.0.4+  Tag:      0.1.0 +Flag withCairo+  Description: Enable non-terminal histogram output, requiring Cairo.+  Default:     False+ Library   HS-Source-Dirs:  src   Ghc-Options:     -O2 -Wall-  Exposed-modules: Statistics.Test.ApproxRand-  Build-Depends:   base >= 4 && < 5, vector == 0.9.*,+  Exposed-modules: Statistics.Test.ApproxRand, Statistics.Test.Correction+  Build-Depends:   base >= 4 && < 5, vector >= 0.9 && < 0.11,                    mersenne-random-pure64 == 0.2.0.*,                    monad-mersenne-random == 0.1,-                   mtl == 2.1.*, statistics == 0.10.1.*,+                   mtl == 2.1.*, statistics == 0.10.*,                    transformers == 0.3.*   Executable approx_rand_test   Main-Is:        approx-rand-test.hs-  Other-Modules:  SampleIO+  Other-Modules:  Histogram, SampleIO, TextHistogram   HS-Source-Dirs: utils   Ghc-Options:    -O2 -Wall   Build-Depends:  base >= 4 && < 5, approx-rand-test,-                  conduit == 0.4.*, text == 0.11.2.*,-                  vector == 0.9.*,+                  conduit >= 0.4 && < 0.6, text == 0.11.2.*,+                  containers == 0.4.*, vector >= 0.9 && < 0.11,                   mersenne-random-pure64 == 0.2.0.*,                   monad-mersenne-random == 0.1,-                  resourcet == 0.3.*,-                  statistics == 0.10.1.*+                  statistics == 0.10.*+  if flag(withCairo)+    Hs-Source-Dirs: utils/cairo+    Build-Depends:  Chart >= 0.15 && < 0.17, colour == 2.3.*,+                    data-accessor == 0.2.2.*, filepath >= 1.3 && < 1.5,+                    vector-algorithms == 0.5.*+  else+    Hs-Source-Dirs: utils/nocairo   Executable approx_rand_test_paired   Main-Is:        approx-rand-test-paired.hs-  Other-Modules:  SampleIO+  Other-Modules:  Histogram, SampleIO, TextHistogram   HS-Source-Dirs: utils   Ghc-Options:    -O2 -Wall   Build-Depends:  base >= 4 && < 5, approx-rand-test,-                  conduit == 0.4.*, text == 0.11.2.*,-                  vector == 0.9.*,+                  conduit >= 0.4 && < 0.6, text == 0.11.2.*,+                  containers == 0.4.*, vector >= 0.9 && < 0.11,                   mersenne-random-pure64 == 0.2.0.*,                   monad-mersenne-random == 0.1,-                  mtl == 2.1.*, resourcet == 0.3.*,-                  statistics == 0.10.1.*+                  mtl == 2.1.*, statistics == 0.10.*+  if flag(withCairo)+    Hs-Source-Dirs: utils/cairo+    Build-Depends:  Chart >= 0.15 && < 0.17, colour == 2.3.*,+                    data-accessor == 0.2.2.*, filepath >= 1.3 && < 1.5,+                    vector-algorithms == 0.5.*++  else+    Hs-Source-Dirs: utils/nocairo  Test-Suite tests   Type:           exitcode-stdio-1.0

@@ -21,15 +21,17 @@   -- $examples    -- * Data types+  TestOptions(..),   TestResult(..),+  Significance(..),   RandWithError,    -- * Approximate randomization tests   approxRandTest,-  approxRandScores,+  approxRandStats,    approxRandPairTest,-  approxRandPairScores,+  approxRandPairStats,    -- * Test statistics   TestStatistic,@@ -45,7 +47,6 @@ import           Control.Monad.Mersenne.Random (R(..), Rand(..), getBool) import           Control.Monad.ST (runST) import           Control.Monad.Trans.Class (lift)-import           Data.List (foldl') import           Data.Vector.Generic ((++)) import qualified Data.Vector.Generic as VG import qualified Data.Vector.Generic.Mutable as GM@@ -97,43 +98,58 @@ -- difference as the test statistic, by running 'approxRandTest' in the 'Rand' -- monad (at the /p = 0.01/ level): ----- > evalRandom (approxRandTest TwoSided meanDifference 10000 0.01 s1 s2) prng+-- > evalRandom (approxRandTest TwoTailed meanDifference 10000 0.01 s1 s2) prng ----- It is also possible to obtain the test scores of the shuffled samples--- directly (e.g. to inspect the distribution of test scores) using the--- 'approxRandScores'/'approxRandPiarScores' functions:+-- It is also possible to obtain the test statistics of the shuffled samples+-- directly (e.g. to inspect the distribution of test statistics) using the+-- 'approxRandStats'/'approxRandPiarStats' functions: ----- > evalRandom (approxRandScores meanDifference 10000 0.01 s1 s2) prng+-- > evalRandom (approxRandStats meanDifference 10000 0.01 s1 s2) prng  -- | Computations with random numbers that can fail. type RandWithError a = ErrorT String Rand a +-- | Options for randomization tests+--+data TestOptions = TestOptions {+    toTestType      :: TestType,      -- ^ Type of test ('OneTailed' or 'TwoTailed')+    toTestStatistic :: TestStatistic, -- ^ Test statistic+    toIterations    :: Int,           -- ^ Number of shuffled samples to create+    toPValue        :: Double         -- ^ he p-value at which to test (e.g. 0.05)+}+ -- | -- The result of hypothesis testing.-data TestResult =-    Significant Double     -- ^ The null hypothesis should be rejected-  | NotSignificant Double  -- ^ Data is compatible with the null hypothesis+data TestResult = TestResult {+    trSignificance     :: Significance, -- ^ Significance+    trStat             :: Double,       -- ^ Test statistic for the samples+    trRandomizedStats  :: Sample        -- ^ Test statistics for the+                                        --   randomized samples+  } deriving (Eq, Ord, Show)++-- |+-- Significance.+data Significance =+    Significant    Double -- ^ The null hypothesis should be rejected+  | NotSignificant Double -- ^ Data is compatible with the null hypothesis   deriving (Eq, Ord, Show)  -- | -- Apply a pair-wise approximate randomization test. ----- In pair-wise approximate randomization tests the scores at a given+-- In pair-wise approximate randomization tests the data points at a given -- index are swapped between samples with a probability of 0.5. Since -- swapping is pairwise, the samples should have the same length. approxRandPairTest ::-     TestType                 -- ^ Type of test ('OneTailed' or 'TwoTailed')-  -> TestStatistic            -- ^ Test statistic-  -> Int                      -- ^ Number of shuffled samples to create-  -> Double                   -- ^ The p-value at which to test (e.g. 0.05)+     TestOptions              -- ^ Options for the test   -> Sample                   -- ^ First sample   -> Sample                   -- ^ Second sample   -> RandWithError TestResult -- ^ The test result-approxRandPairTest testType stat n pTest s1 s2 =-  (significance testType pTest n . countExtremes tOrig) `liftM`-    approxRandPairScores stat n s1 s2-  where-    tOrig = stat s1 s2+approxRandPairTest (TestOptions testType stat n pTest) s1 s2 = do+  stats <- approxRandPairStats stat n s1 s2+  let tOrig = stat s1 s2+  let sig = significance testType pTest n $ countExtremes tOrig $ stats+  return $ TestResult sig tOrig stats  -- | -- Apply an approximate randomization test.@@ -143,26 +159,23 @@ -- the original samples and the shuffled samples, to detect whether the -- difference of the samples is extreme or not. approxRandTest ::-     TestType        -- ^ Type of test ('OneTailed' or 'TwoTailed')-  -> TestStatistic   -- ^ Test statistic-  -> Int             -- ^ Number of shuffled sample to create-  -> Double          -- ^ The p-value at which to test (e.g. 0.05)+     TestOptions     -- ^ Options for the test   -> Sample          -- ^ First sample   -> Sample          -- ^ Second sample   -> Rand TestResult -- ^ The test result-approxRandTest testType stat n pTest s1 s2 =-  (significance testType pTest n . countExtremes tOrig) `liftM`-    approxRandScores stat n s1 s2-  where-    tOrig = stat s1 s2+approxRandTest (TestOptions testType stat n pTest) s1 s2 = do+  stats <- approxRandStats stat n s1 s2+  let tOrig = stat s1 s2+  let sig = significance testType pTest n $ countExtremes tOrig stats+  return $ TestResult sig tOrig stats  -- | Determine the significance. significance ::-     TestType   -- ^ Type of test ('OneTailed' or 'TwoTailed')-  -> Double     -- ^ The p-value at which to test (e.g. 0.05)-  -> Int        -- ^ Number of sample shuffles-  -> (Int, Int) -- ^ Extreme score counts-  -> TestResult -- ^ The test result+     TestType     -- ^ Type of test ('OneTailed' or 'TwoTailed')+  -> Double       -- ^ The p-value at which to test (e.g. 0.05)+  -> Int          -- ^ Number of sample shuffles+  -> (Int, Int)   -- ^ Extreme statistic counts+  -> Significance -- ^ The test result significance TwoTailed pTest n =   significant (pTest / 2) . pValue n . uncurry min significance OneTailed pTest n =@@ -170,15 +183,15 @@  -- | Wrap a p-value in a 'TestResult'. significant ::-     Double     -- ^ The p-value at which to test-  -> Double     -- ^ The p-value-  -> TestResult -- ^ The test result+     Double       -- ^ The p-value at which to test+  -> Double       -- ^ The p-value+  -> Significance -- ^ The test result significant pTail p =   if p < pTail then Significant p else NotSignificant p  -- | Calculate a p-value pValue ::-     Int    -- ^ Number of extreme scores+     Int    -- ^ Number of extreme outcomes   -> Int    -- ^ Number of shuffles   -> Double -- ^ The p-value pValue n r = (fromIntegral r + 1) / (fromIntegral n + 1)@@ -190,14 +203,14 @@ -- count value smaller than or equal to that value. Since we do not know -- the tail (yet), we count both. ----- Note: we can determine the tail by (1) averaging the test scores of the--- randomized samples, or (2) taking the smaller of the two counts.+-- Note: we can determine the tail by (1) averaging the test statistics of+-- the randomized samples, or (2) taking the smaller of the two counts. countExtremes ::      Double     -- ^ Test statistic value of the original samples-  -> [Double]   -- ^ Test statistic values of the randomized samples.+  -> Sample     -- ^ Test statistic values of the randomized samples.   -> (Int, Int) -- ^ Count of left- and right-tail extremes. countExtremes tOrig =-  foldl' count (0, 0)+  VG.foldl' count (0, 0)   where     count (left, right) tPerm =       let !newLeft = if tPerm <= tOrig then succ left else left in@@ -206,34 +219,35 @@  -- | -- Generate a given number of pairwise shuffled samples, and calculate--- the test score for each shuffle.+-- the test statistic for each shuffle. ----- Since the scores at a given index are swapped (with a probability of+-- Since the data points at a given index are swapped (with a probability of -- 0.5), the samples should have the same length.-approxRandPairScores ::+approxRandPairStats ::      TestStatistic          -- ^ Test statistic   -> Int                    -- ^ Number of shuffled samples to create   -> Sample                 -- ^ First sample   -> Sample                 -- ^ Second sample-  -> RandWithError [Double] -- ^ The scores of each shuffle-approxRandPairScores stat n s1 s2 = do+  -> RandWithError Sample   -- ^ The statistics of the shuffles+approxRandPairStats stat n s1 s2 = do   when (VG.length s1 /= VG.length s2) $-    throwError "Cannot calculate pairwise scores: samples have different sizes"-  lift $ replicateM n $ uncurry stat `liftM` shuffleVectorsPairwise s1 s2+    throwError "Cannot calculate pairwise statistic: samples have different sizes"+  lift $ liftM VG.fromList $ replicateM n $+    uncurry stat `liftM` shuffleVectorsPairwise s1 s2  -- | -- Generate a given number of shuffled samples, and calculate the test--- score for each shuffle.+-- statistic for each shuffle. -- -- This function does not require the samples to have an equal length.-approxRandScores ::+approxRandStats ::      TestStatistic -- ^ Test statistic   -> Int           -- ^ Number of shuffled samples to create   -> Sample        -- ^ First sample   -> Sample        -- ^ Second sample-  -> Rand [Double] -- ^ The scores of each shuffle-approxRandScores stat n s1 s2 =-  replicateM n $ uncurry stat `liftM` shuffleVectors s1 s2+  -> Rand Sample   -- ^ The statistics of the shuffles+approxRandStats stat n s1 s2 =+  liftM VG.fromList $ replicateM n $ uncurry stat `liftM` shuffleVectors s1 s2  -- | Pair-wise shuffle of two vectors. shuffleVectorsPairwise :: (VG.Vector v a, VG.Vector v Bool) =>

@@ -0,0 +1,42 @@+-- |+-- Copyright  : (c) 2012 Daniël de Kok+-- License    : BSD3+--+-- Maintainer : Daniël de Kok <me@danieldk.eu>+-- Stability  : experimental+--+-- The more tests that are applied on a dataset, the more likely it becomes+-- that a Type I error occurs (rejecting the null hypothesis when it is+-- true). This module provides significance level corrections for repeated+-- tests.+--+-- For more information, see:+--+--  * /Bonferroni and Šidák corrections for multiple comparisons, H. Abdi, 2007, in: N.J. Salkind (ed.), Encyclopedia of Measurement and Statistics, Thousand Oaks, CA: Sage/++module Statistics.Test.Correction (+  sidak,+  bonferroni+) where++-- |+-- Bonferroni correction: α* = α / n, where α is the intended+-- significance level over all tests, n the number of tests, and α* the+-- corrected significance level.+bonferroni :: (Fractional f, Integral i) =>+     f -- ^ Significance level+  -> i -- ^ Number of comparisons+  -> f -- ^ The corrected significance level+bonferroni level n =+  level / fromIntegral n++-- |+-- Sidak correction: α* = 1 - (1 - α)^(1 / n), where α is the intended+-- significance level over all tests, n the number of tests, and α* the+-- corrected significance level.+sidak :: (Floating f, Integral i) =>+     f -- ^ Significance level+  -> i -- ^ Number of comparisons+  -> f -- ^ The corrected significance level+sidak level n =+  1 - (1 - level) ** (1 / (fromIntegral n))

@@ -0,0 +1,33 @@+module Histogram (+  histogram+) where++import qualified Data.Set as Set+import qualified Data.Vector.Generic as VG+import           Data.Vector.Unboxed ((!))+import qualified Statistics.Sample.Histogram as SSH+import           Statistics.Test.ApproxRand+import           Statistics.Types (Sample)++histogram ::+     Int+  -> TestResult+  -> Either String [(Double, Int)]+histogram prefBins (TestResult _ _ randomizedStats) =+  if bins < 2 then+    Left "Cannot make a histogram, because the data cannot be divided in more than one bin."+  else+    Right $ VG.toList $ VG.zip ticks numSamples+  where+    bins = numBins prefBins randomizedStats+    (lowerBounds, numSamples) =+      SSH.histogram prefBins randomizedStats++    -- Compute the middle values of the bins.+    bucketSize = (lowerBounds ! 1) - (lowerBounds ! 0)+    bucketHalf = bucketSize / 2+    ticks = VG.map ((+) bucketHalf) lowerBounds++numBins :: Int -> Sample -> Int+numBins pref =+  min pref . Set.size . VG.foldl' (flip Set.insert) Set.empty

@@ -0,0 +1,30 @@+module TextHistogram (+  printHistogram,+) where++import           Control.Monad (forM_, when)+import           Statistics.Test.ApproxRand+import           System.IO (hPutStrLn, stderr)+import           Text.Printf (printf)++import           Histogram++printHistogram :: Int -> TestResult -> IO ()+printHistogram bins result@(TestResult _ score _) =+  case histogram bins result of+    Left err   -> hPutStrLn stderr err+    Right hist ->+      let bucketSize = (fst $ hist !! 1) - (fst $ head hist)+          bucketHalf = bucketSize / 2+          charsPerDot = (maximum $ map snd hist) `div` 50 in+            forM_ hist $ \(label, freq) -> do+              let blocks = freq `div` charsPerDot+              when (blocks > 0) $ do+                let lower = label - bucketHalf+                let barChar = if score >= lower && score < lower + bucketSize then+                            '✣'+                          else+                            '█'+                putStr $ printf "%12.3e | " $ lower + bucketHalf+                putStrLn $ replicate blocks barChar+

@@ -14,6 +14,7 @@ import           Control.Monad (liftM, when) import           Control.Monad.Error (runErrorT) import           Control.Monad.Mersenne.Random (evalRandom)+import qualified Data.Vector.Generic as VG import qualified Data.Vector.Unboxed as V import           Data.Word (Word64) import           Statistics.Test.ApproxRand@@ -25,7 +26,9 @@ import           System.Random.Mersenne.Pure64 (PureMT, newPureMT, pureMT) import           Text.Printf (printf) +import           CairoHistogram import           SampleIO+import           TextHistogram  main :: IO () main = do@@ -43,8 +46,8 @@     Just seed -> return $ pureMT seed     Nothing   -> newPureMT -  if optPrintScores opts then-    printScores opts stat prng v1 v2+  if optPrintStats opts then+    printStats opts stat prng v1 v2   else     applyTest opts stat prng v1 v2 @@ -55,9 +58,10 @@   putStrLn $ printf "Iterations: %d" $ optIterations opts   putStrLn $ printf "Sample size: %d" $ V.length v1 -  -- Calculate test statistic for original score sets.-  let tOrig = stat v1 v2-  putStrLn $ printf "Test statistic: %f" tOrig+  let testOptions = TestOptions (optTestType opts)+                                stat+                                (optIterations opts)+                                (optSigP opts)    let testType = optTestType opts @@ -71,48 +75,81 @@   putStrLn $ printf "Test significance: %f" pTest   putStrLn $ printf "Tail significance: %f" pTail -  -- Approximate randomization testing.-  let test = runErrorT $ approxRandPairTest testType stat (optIterations opts) pTest v1 v2+  let test = runErrorT $ approxRandPairTest testOptions v1 v2   let result = evalRandom test prng+   case result of-    Left  err                -> putStrLn err-    Right (Significant    p) -> putStrLn $ printf "Significant: %f" p-    Right (NotSignificant p) -> putStrLn $ printf "Not significant: %f" p+    Left  err -> putStrLn err+    Right r   -> do+                   printResult r+                   when (optPrintHistogram opts) $ do+                     putStrLn ""+                     printHistogram 21 r+                   case (optWriteHistogram opts) of+                     Just fn ->+                       writeHistogram testOptions 31 r fn+                     Nothing ->+                       return () -printScores :: Options -> TestStatistic -> PureMT -> Sample ->+printResult :: TestResult -> IO ()+printResult result = do+  -- Print test statistic for the samples.+  putStrLn $ printf "Test statistic: %f" $ trStat result++  -- Print significance+  case trSignificance result of+    Significant    p -> putStrLn $ printf "Significant: %f" p+    NotSignificant p -> putStrLn $ printf "Not significant: %f" p++printStats :: Options -> TestStatistic -> PureMT -> Sample ->   Sample -> IO ()-printScores opts stat prng v1 v2 = do-  let test = runErrorT $ approxRandPairScores stat (optIterations opts) v1 v2+printStats opts stat prng v1 v2 = do+  let test = runErrorT $ approxRandPairStats stat (optIterations opts) v1 v2   case evalRandom test prng of     Left err     -> putStrLn err-    Right scores -> mapM_ (putStrLn . printf "%f") scores+    Right scores -> VG.mapM_ (putStrLn . printf "%f") scores  data Options = Options {-  optColumn        :: Int,-  optIterations    :: Int,-  optPRNGSeed      :: Maybe Word64,-  optPrintScores   :: Bool,-  optSigP          :: Double,-  optTestStatistic :: TestStatistic,-  optTestType      :: TestType+  optColumn         :: Int,+  optIterations     :: Int,+  optPRNGSeed       :: Maybe Word64,+  optPrintHistogram :: Bool,+  optPrintStats     :: Bool,+  optSigP           :: Double,+  optTestStatistic  :: TestStatistic,+  optTestType       :: TestType,+  optWriteHistogram :: Maybe String }  defaultOptions :: Options defaultOptions = Options {-  optColumn        = 1,-  optIterations    = 10000,-  optPRNGSeed      = Nothing,-  optPrintScores   = False,-  optSigP          = 0.01,-  optTestStatistic = differenceMean,-  optTestType      = TwoTailed+  optColumn         = 1,+  optIterations     = 10000,+  optPRNGSeed       = Nothing,+  optPrintHistogram = False,+  optPrintStats     = False,+  optSigP           = 0.01,+  optTestStatistic  = differenceMean,+  optTestType       = TwoTailed,+  optWriteHistogram = Nothing }  options :: [OptDescr (Options -> Options)] options =+  if hasCairoHistograms then+    cairoHistogramOption : mandatoryOptions+  else+    mandatoryOptions++-- Options that are always available, regardless of the compilation flags.+mandatoryOptions :: [OptDescr (Options -> Options)]+mandatoryOptions =   [ Option ['c'] ["column"]       (ReqArg (\arg opt -> opt { optColumn = read arg }) "NUMBER")       "column number (starting at 1)",+    Option ['h']    ["print-histogram"]+      (NoArg (\opt -> opt { optPrintHistogram = True }))+      "print a histogram of randomized sample statistics",     Option ['i'] ["iterations"]       (ReqArg (\arg opt -> opt { optIterations = read arg }) "NUMBER")       "number of iterations",@@ -122,9 +159,9 @@     Option ['p'] []       (ReqArg (\arg opt -> opt {optSigP = read arg }) "NUMBER")       "significant p-value",-    Option []    ["print-scores"]-      (NoArg (\opt -> opt { optPrintScores = True }))-      "output scores of permuted vectors",+    Option []    ["print-statistics"]+      (NoArg (\opt -> opt { optPrintStats = True }))+      "output test statistics of permuted vectors",     Option ['s'] ["seed"]       (ReqArg (\arg opt -> opt { optPRNGSeed = Just $ read arg}) "NUMBER")       "pseudorandom number generator seed",@@ -132,6 +169,12 @@       (ReqArg (\arg opt -> opt { optTestStatistic = parseStatistic arg}) "NAME")       "test statistic (mean_diff, var_ratio)"   ]++cairoHistogramOption :: OptDescr (Options -> Options)+cairoHistogramOption =+    Option ['w'] ["write-histogram"]+      (ReqArg (\arg opt -> opt { optWriteHistogram = Just arg}) "FILENAME")+      "write a histogram (supported file extensions: pdf, png, ps, and svg)"  getOptions :: IO (Options, [String]) getOptions = do

@@ -13,6 +13,7 @@  import           Control.Monad (liftM, when) import           Control.Monad.Mersenne.Random (evalRandom)+import qualified Data.Vector.Generic as VG import qualified Data.Vector.Unboxed as V import           Data.Word (Word64) import           Statistics.Test.ApproxRand@@ -24,7 +25,9 @@ import           System.Random.Mersenne.Pure64 (PureMT, newPureMT, pureMT) import           Text.Printf (printf) +import           CairoHistogram import           SampleIO+import           TextHistogram  main :: IO () main = do@@ -42,8 +45,8 @@     Just seed -> return $ pureMT seed     Nothing   -> newPureMT -  if optPrintScores opts then-    printScores opts stat prng v1 v2+  if optPrintStats opts then+    printStats opts stat prng v1 v2   else     applyTest opts stat prng v1 v2 @@ -54,9 +57,10 @@   putStrLn $ printf "Iterations: %d" $ optIterations opts   putStrLn $ printf "Sample sizes: %d %d" (V.length v1) (V.length v2) -  -- Calculate test statistic for original score sets.-  let tOrig = stat v1 v2-  putStrLn $ printf "Test statistic: %f" tOrig+  let testOptions = TestOptions (optTestType opts)+                                stat+                                (optIterations opts)+                                (optSigP opts)    let testType = optTestType opts @@ -71,44 +75,71 @@   putStrLn $ printf "Tail significance: %f" pTail    -- Approximate randomization testing.-  let test = approxRandTest testType stat (optIterations opts) pTest v1 v2+  let test = approxRandTest testOptions v1 v2   let result = evalRandom test prng-  case result of++  -- Print test statistic for the samples.+  putStrLn $ printf "Test statistic: %f" $ trStat result++  case trSignificance result of     Significant    p -> putStrLn $ printf "Significant: %f" p     NotSignificant p -> putStrLn $ printf "Not significant: %f" p -printScores :: Options -> TestStatistic -> PureMT -> Sample ->+  when (optPrintHistogram opts) $ do+    putStrLn ""+    printHistogram 21 result+  case (optWriteHistogram opts) of+    Just fn ->+      writeHistogram testOptions 31 result fn+    Nothing ->+      return ()++printStats :: Options -> TestStatistic -> PureMT -> Sample ->   Sample -> IO ()-printScores opts stat prng v1 v2 =-  mapM_ (putStrLn . printf "%f") $-    evalRandom (approxRandScores stat (optIterations opts) v1 v2) prng+printStats opts stat prng v1 v2 =+  VG.mapM_ (putStrLn . printf "%f") $+    evalRandom (approxRandStats stat (optIterations opts) v1 v2) prng  data Options = Options {-  optColumn        :: Int,-  optIterations    :: Int,-  optPRNGSeed      :: Maybe Word64,-  optPrintScores   :: Bool,-  optSigP          :: Double,-  optTestStatistic :: TestStatistic,-  optTestType      :: TestType+  optColumn         :: Int,+  optPrintHistogram :: Bool,+  optIterations     :: Int,+  optPRNGSeed       :: Maybe Word64,+  optPrintStats     :: Bool,+  optSigP           :: Double,+  optTestStatistic  :: TestStatistic,+  optTestType       :: TestType,+  optWriteHistogram :: Maybe String }  defaultOptions :: Options defaultOptions = Options {-  optColumn        = 1,-  optIterations    = 10000,-  optPRNGSeed      = Nothing,-  optPrintScores   = False,-  optSigP          = 0.01,-  optTestStatistic = meanDifference,-  optTestType      = TwoTailed+  optColumn         = 1,+  optIterations     = 10000,+  optPRNGSeed       = Nothing,+  optPrintHistogram = False,+  optPrintStats     = False,+  optSigP           = 0.01,+  optTestStatistic  = meanDifference,+  optTestType       = TwoTailed,+  optWriteHistogram = Nothing }  options :: [OptDescr (Options -> Options)]-options =+options = +  if hasCairoHistograms then+    cairoHistogramOption : mandatoryOptions+  else+    mandatoryOptions++mandatoryOptions :: [OptDescr (Options -> Options)]+mandatoryOptions =   [ Option ['c'] ["column"]       (ReqArg (\arg opt -> opt { optColumn = read arg }) "NUMBER")       "column number (starting at 1)",+    Option ['h']    ["print-histogram"]+      (NoArg (\opt -> opt { optPrintHistogram = True }))+      "print a histogram of randomized sample statistics",     Option ['i'] ["iterations"]       (ReqArg (\arg opt -> opt { optIterations = read arg }) "NUMBER")       "number of iterations",@@ -119,8 +150,8 @@       (ReqArg (\arg opt -> opt {optSigP = read arg }) "NUMBER")       "significant p-value",     Option []    ["print-scores"]-      (NoArg (\opt -> opt { optPrintScores = True }))-      "output scores of permuted vectors",+      (NoArg (\opt -> opt { optPrintStats = True }))+      "output statistics of permuted vectors",     Option ['s'] ["seed"]       (ReqArg (\arg opt -> opt { optPRNGSeed = Just $ read arg}) "NUMBER")       "pseudorandom number generator seed",@@ -128,6 +159,12 @@       (ReqArg (\arg opt -> opt { optTestStatistic = parseStatistic arg}) "NAME")       "test statistic (mean_diff, var_ratio)"   ]++cairoHistogramOption :: OptDescr (Options -> Options)+cairoHistogramOption =+    Option ['w'] ["write-histogram"]+      (ReqArg (\arg opt -> opt { optWriteHistogram = Just arg}) "FILENAME")+      "write a histogram (supported file extensions: pdf, png, ps, and svg)"  getOptions :: IO (Options, [String]) getOptions = do

@@ -0,0 +1,104 @@+module CairoHistogram (+  hasCairoHistograms,+  writeHistogram+) where++import           Control.Monad.ST (runST)+import           Data.Accessor ((^=), (^:))+import qualified Data.Colour as Colour+import qualified Data.Colour.Names as ColourNames+import qualified Data.Vector.Algorithms.Intro as VI+import qualified Data.Vector.Generic as VG+import           Data.Vector.Unboxed ((!))+import qualified Graphics.Rendering.Chart as Chart+import           Statistics.Test.ApproxRand+import           Statistics.Test.Types (TestType(..))+import qualified System.FilePath.Posix as FP+import           System.IO (hPutStrLn, stderr)++import           Histogram++hasCairoHistograms :: Bool+hasCairoHistograms = True++writeHistogram :: TestOptions -> Int -> TestResult -> FP.FilePath -> IO ()+writeHistogram testOptions bins result path =+  case histogram bins result of+    Left err -> hPutStrLn stderr err+    Right  h ->+      let r = createHistogram testOptions result h in+        case snd $ FP.splitExtension path of+          ".pdf" -> Chart.renderableToPDFFile r 800 600 path+          ".png" -> do+            _ <- Chart.renderableToPNGFile r 800 600 path+            return ()+          ".ps"  -> Chart.renderableToPSFile  r 800 600 path+          ".svg" -> Chart.renderableToSVGFile r 800 600 path+          _      -> hPutStrLn stderr "Unknown output format!"++-- Creates a histogram. The histogram is stacked, but the second bar+-- is always empty, except for the bin of the original statistic (if any).+-- There, the first bar is empty and the second bar has the frequency.+-- Yes, this is cheating ;).+createHistogram :: TestOptions -> TestResult -> [(Double, Int)] ->+  Chart.Renderable ()+createHistogram testOptions result his =+  Chart.toRenderable layout+  where+    layout =+        Chart.layout1_background  ^= Chart.solidFillStyle opaqueWhite+      $ Chart.layout1_left_axis   ^: Chart.laxis_override ^= Chart.axisTicksHide+      $ Chart.layout1_right_axis  ^: Chart.laxis_title    ^= "Frequency"+      $ Chart.layout1_bottom_axis ^: Chart.laxis_title    ^= "Statistic"+      $ Chart.layout1_plots       ^= [ Right (Chart.plotBars randomizationBars),+                                       Right statisticLine, Right sigLines ]+      $ Chart.setLayout1Foreground   opaqueBlack+      $ Chart.defaultLayout1+    randomizationBars =+        Chart.plot_bars_style       ^= Chart.BarsStacked+      $ Chart.plot_bars_spacing     ^= Chart.BarsFixGap 6 2+      -- $ Chart.plot_bars_spacing     ^= Chart.BarsFixGap 0 0+      $ Chart.plot_bars_item_styles ^= [+          (Chart.solidFillStyle $ opaqueGreen, Nothing) ]+      $ Chart.plot_bars_values      ^= map (\(b, f) -> (b, [f])) his+      $ Chart.defaultPlotBars+    statisticLine =+      Chart.vlinePlot "Statistic for samples" (Chart.solidLine 2 (opaqueRed)) $ trStat result+    sigLines      =+      vlinesPlot "Significance" (Chart.dashedLine 2 [8, 4] opaqueBlack) $+        sigBounds testOptions result++-- Plot vertical lines, adapted from Chart.vlinePlot.+vlinesPlot :: String -> Chart.CairoLineStyle -> [a] -> Chart.Plot a b+vlinesPlot t ls xs = Chart.toPlot Chart.defaultPlotLines {+    Chart.plot_lines_title_        = t,+    Chart.plot_lines_style_        = ls,+    Chart.plot_lines_limit_values_ =+      [[(Chart.LValue v, Chart.LMin),(Chart.LValue v, Chart.LMax)] | v <- xs]+    }++-- Calculate the bounds of significance.+sigBounds :: TestOptions -> TestResult -> [Double]+sigBounds (TestOptions testType _ n pTest) (TestResult _ _ stats) =+  case testType of+    TwoTailed -> [sorted ! (nExtreme - 1), sorted ! (n - nExtreme)]+    OneTailed -> [sorted ! (n - nExtreme)]+  where+    sorted           = sortVector stats+    nExtreme         = floor $ (pVal testType pTest) * (fromIntegral n + 1) - 1+    pVal OneTailed p = p+    pVal TwoTailed p = p / 2+    -- XXX: Fix extreme cases: p-value of 0, small n.++sortVector :: (Ord a, VG.Vector v a) => v a -> v a+sortVector v = runST $ do+  s <- VG.thaw v+  VI.sort s+  VG.freeze s++-- Convenience...+opaqueBlack, opaqueGreen, opaqueRed, opaqueWhite :: Colour.AlphaColour Double+opaqueBlack = Colour.opaque ColourNames.black+opaqueGreen = Colour.opaque ColourNames.green+opaqueRed   = Colour.opaque ColourNames.red+opaqueWhite = Colour.opaque ColourNames.white

@@ -0,0 +1,15 @@+module CairoHistogram (+  hasCairoHistograms,+  writeHistogram+) where++import           System.IO (hPutStrLn, stderr)++hasCairoHistograms :: Bool+hasCairoHistograms = False++-- Stub function, for when we did not compile with Cairo support.+writeHistogram :: a -> Int -> b -> c -> IO ()+writeHistogram _ _ _ _ = do+  hPutStrLn stderr "This program was built without Cairo support, cannot make a histogram."+