diff --git a/approx-rand-test.cabal b/approx-rand-test.cabal --- a/approx-rand-test.cabal +++ b/approx-rand-test.cabal @@ -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 diff --git a/src/Statistics/Test/ApproxRand.hs b/src/Statistics/Test/ApproxRand.hs --- a/src/Statistics/Test/ApproxRand.hs +++ b/src/Statistics/Test/ApproxRand.hs @@ -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) => diff --git a/src/Statistics/Test/Correction.hs b/src/Statistics/Test/Correction.hs new file mode 100644 --- /dev/null +++ b/src/Statistics/Test/Correction.hs @@ -0,0 +1,42 @@ +-- | +-- Copyright : (c) 2012 Daniël de Kok +-- License : BSD3 +-- +-- Maintainer : Daniël de Kok +-- 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)) diff --git a/utils/Histogram.hs b/utils/Histogram.hs new file mode 100644 --- /dev/null +++ b/utils/Histogram.hs @@ -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 diff --git a/utils/TextHistogram.hs b/utils/TextHistogram.hs new file mode 100644 --- /dev/null +++ b/utils/TextHistogram.hs @@ -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 + diff --git a/utils/approx-rand-test-paired.hs b/utils/approx-rand-test-paired.hs --- a/utils/approx-rand-test-paired.hs +++ b/utils/approx-rand-test-paired.hs @@ -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 diff --git a/utils/approx-rand-test.hs b/utils/approx-rand-test.hs --- a/utils/approx-rand-test.hs +++ b/utils/approx-rand-test.hs @@ -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 diff --git a/utils/cairo/CairoHistogram.hs b/utils/cairo/CairoHistogram.hs new file mode 100644 --- /dev/null +++ b/utils/cairo/CairoHistogram.hs @@ -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 diff --git a/utils/nocairo/CairoHistogram.hs b/utils/nocairo/CairoHistogram.hs new file mode 100644 --- /dev/null +++ b/utils/nocairo/CairoHistogram.hs @@ -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." +