diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,202 @@ + + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,4 @@ +#!/usr/bin/env runhaskell + +> import Distribution.Simple +> main = defaultMain diff --git a/approx-rand-test.cabal b/approx-rand-test.cabal new file mode 100644 --- /dev/null +++ b/approx-rand-test.cabal @@ -0,0 +1,71 @@ +Name: approx-rand-test +Version: 0.0.3 +License: OtherLicense +License-file: LICENSE +Copyright: Copyright 2012 Daniël de Kok +Author: Daniël de Kok +Maintainer: Daniël de Kok +Homepage: http://github.com/danieldk/approx-rand-test +Category: Statistics +Synopsis: Approximate randomization test +Description: Utility to perform approximate randomization tests. +Cabal-Version: >= 1.8 +Build-Type: Simple + +Source-Repository head + Type: git + Location: git://github.com/danieldk/approx-rand-test.git + +Source-Repository this + Type: git + Location: git://github.com/danieldk/approx-rand-test.git + Tag: 0.0.3 + +Library + HS-Source-Dirs: src + Ghc-Options: -O2 -Wall + Exposed-modules: Statistics.Test.ApproxRand + Build-Depends: base >= 4 && < 5, vector == 0.9.*, + mersenne-random-pure64 == 0.2.0.*, + monad-mersenne-random == 0.1, + mtl == 2.0.1.*, statistics == 0.10.1.*, + transformers == 0.2.2.* + + +Executable approx_rand_test + Main-Is: approx-rand-test.hs + HS-Source-Dirs: utils + Ghc-Options: -O2 -Wall + Build-Depends: base >= 4 && < 5, approx-rand-test, + conduit == 0.4.*, text == 0.11.1.*, + vector == 0.9.*, + mersenne-random-pure64 == 0.2.0.*, + monad-mersenne-random == 0.1, + resourcet == 0.3.*, + statistics == 0.10.1.* + + +Executable approx_rand_test_paired + Main-Is: approx-rand-test-paired.hs + HS-Source-Dirs: utils + Ghc-Options: -O2 -Wall + Build-Depends: base >= 4 && < 5, approx-rand-test, + conduit == 0.4.*, text == 0.11.1.*, + vector == 0.9.*, + mersenne-random-pure64 == 0.2.0.*, + monad-mersenne-random == 0.1, + mtl == 2.0.1.*, resourcet == 0.3.*, + statistics == 0.10.1.* + +Test-Suite tests + Type: exitcode-stdio-1.0 + Hs-Source-Dirs: tests + Main-Is: tests.hs + ghc-options: -Wall + Build-Depends: base >= 4 && < 5, vector == 0.9.*, + approx-rand-test, + mersenne-random-pure64 == 0.2.0.*, + monad-mersenne-random == 0.1, + HUnit == 1.2.4.*, + test-framework == 0.6.*, + test-framework-hunit == 0.2.* diff --git a/src/Statistics/Test/ApproxRand.hs b/src/Statistics/Test/ApproxRand.hs new file mode 100644 --- /dev/null +++ b/src/Statistics/Test/ApproxRand.hs @@ -0,0 +1,342 @@ +-- | +-- Copyright : (c) 2012 Daniël de Kok +-- License : BSD3 +-- +-- Maintainer : Daniël de Kok +-- Stability : experimental +-- +-- This module provides functionality to perform approximate randomization +-- tests (Noreen, 1989). + + +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE DoAndIfThenElse #-} +{-# LANGUAGE FlexibleContexts #-} + +module Statistics.Test.ApproxRand ( + -- * Description + -- $description + + -- * Examples + -- $examples + + -- * Data types + TestResult(..), + RandWithError, + + -- * Approximate randomization tests + approxRandTest, + approxRandScores, + + approxRandPairTest, + approxRandPairScores, + + -- * Test statistics + TestStatistic, + differenceMean, + meanDifference, + varianceRatio +) where + +import Prelude hiding ((++)) +import Control.Monad (liftM, replicateM, when) +import Control.Monad.Error (ErrorT) +import Control.Monad.Error.Class (throwError) +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 +import Data.Word (Word) +import Statistics.Sample (variance) +import Statistics.Test.Types (TestType(..)) +import Statistics.Types +import System.Random.Mersenne.Pure64 (PureMT, randomInt, randomWord) + + +-- $description +-- +-- Approximate randomization tests rely on a simple premise: given a test +-- statistic, if the null-hypothesis (the samples do not differ) is true, +-- we can randomly swap values between samples without an (extreme) impact +-- on the test statistic. Otherwise, the null-hypothesis must be rejected. +-- +-- The test works by generating a given number of sample shuffles and computing +-- the test statistic for each shuffle. If /r/ is the number of shuffled +-- samples where the test statistic is at least as high as the test statistic +-- applied on the original samples; and /N/ the number of shuffles, then +-- the null-hypothesis is rejected iff /(r + 1):(N + 1) < p-value/ (for +-- one-sided tests). +-- +-- Two kinds of test are supported: +-- +-- * /Paired sample/ ('approxRandPairTest'): values from samples are shuffled +-- pair-wise. This requires the samples to have an equal length. +-- +-- * /Unpaired sample/ ('approxRandTest'): values from samples are shuffled +-- among both samples. Consequently the i-th element of one sample does not +-- bear a relationship with the i-th element of the other sample. The +-- shuffled samples retain the sizes of the original samples. +-- +-- Both tests can be performed as a one-tailed or two-tailed test. + +-- $examples +-- Both unpaired and paired sample tests use the 'Rand' monad to obtain +-- random numbers. We can obtain a pseudo-random number generator that +-- is seeded using the system clock using the +-- 'System.Random.Mersenne.Pure64.newPureMT' function (please refer to +-- the documentation of 'System.Random.Mersenne.Pure64' for more +-- information): +-- +-- > prng <- newPureMT +-- +-- Suppose that we have the samples 's1' and 's2'. We could now perform +-- a Two-Tailed randomization test with 10,000 shuffles and the mean +-- 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 +-- +-- 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: +-- +-- > evalRandom (approxRandScores meanDifference 10000 0.01 s1 s2) prng + +-- | Computations with random numbers that can fail. +type RandWithError a = ErrorT String Rand a + +-- | +-- The result of hypothesis testing. +data TestResult = + 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 +-- 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) + -> 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 + +-- | +-- Apply an approximate randomization test. +-- +-- In approximate randomization tests, the values of two samples are +-- shuffled among those samples. A test statistic is calculated for +-- 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) + -> 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 + +-- | 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 +significance TwoTailed pTest n = + significant (pTest / 2) . pValue n . uncurry min +significance OneTailed pTest n = + significant pTest . pValue n . snd + +-- | Wrap a p-value in a 'TestResult'. +significant :: + Double -- ^ The p-value at which to test + -> Double -- ^ The p-value + -> TestResult -- ^ 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 shuffles + -> Double -- ^ The p-value +pValue n r = (fromIntegral r + 1) / (fromIntegral n + 1) + +-- | +-- Count extreme test statistic values. If the test statistic value of the +-- original sample is in the right tail, we want to count values equal to +-- or larger than that value. If the value is in the left tail, we want to +-- 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. +countExtremes :: + Double -- ^ Test statistic value of the original samples + -> [Double] -- ^ Test statistic values of the randomized samples. + -> (Int, Int) -- ^ Count of left- and right-tail extremes. +countExtremes tOrig = + foldl' count (0, 0) + where + count (left, right) tPerm = + let !newLeft = if tPerm <= tOrig then succ left else left in + let !newRight = if tPerm >= tOrig then succ right else right in + (newLeft, newRight) + +-- | +-- Generate a given number of pairwise shuffled samples, and calculate +-- the test score for each shuffle. +-- +-- Since the scores at a given index are swapped (with a probability of +-- 0.5), the samples should have the same length. +approxRandPairScores :: + 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 + 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 + +-- | +-- Generate a given number of shuffled samples, and calculate the test +-- score for each shuffle. +-- +-- This function does not require the samples to have an equal length. +approxRandScores :: + 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 + +-- | Pair-wise shuffle of two vectors. +shuffleVectorsPairwise :: (VG.Vector v a, VG.Vector v Bool) => + v a -> v a -> Rand (v a, v a) +shuffleVectorsPairwise vec1 vec2 = do + randomVec <- randomVector (VG.length vec1) + let pv1 = VG.zipWith3 permute vec1 vec2 randomVec + let pv2 = VG.zipWith3 permute vec2 vec1 randomVec + return (pv1, pv2) + where + permute val1 val2 coin = + if coin then val1 else val2 + +randomVector :: (VG.Vector v Bool) => Int -> Rand (v Bool) +randomVector len = + VG.replicateM len getBool + +-- Shuffle values amongst two vectors, keeping the original vector lengths. +shuffleVectors :: VG.Vector v a => v a -> v a -> Rand (v a, v a) +shuffleVectors v1 v2 = do + shuffledVectors <- shuffleVector $ v1 ++ v2 + return (VG.slice 0 (VG.length v1) shuffledVectors, + VG.slice (VG.length v1) (VG.length v2) shuffledVectors) + +-- Fisher-Yates shuffle in the Rand monad +shuffleVector :: VG.Vector v a => v a -> Rand (v a) +shuffleVector v = + Rand $ \s -> case shuffleVector' s v of (sv, s') -> R sv s' + +-- Fisher-Yates shuffle +shuffleVector' :: VG.Vector v a => PureMT -> v a -> (v a, PureMT) +shuffleVector' gen v = runST $ do + let maxIdx = VG.length v - 1 + vm <- VG.thaw v + gen' <- swaps vm 0 maxIdx gen + vmf <- VG.unsafeFreeze vm + return (vmf, gen') + where + swaps vm idx maxIdx gen' + | idx < maxIdx = do + let (newIdx, gen'') = randomIntR gen' (idx, maxIdx) + GM.unsafeSwap vm idx newIdx + swaps vm (idx + 1) maxIdx gen'' + | otherwise = return gen' + +-- | +-- A test stastic calculates the difference between two samples. See +-- 'meanDifference' and 'varianceRatio' for examples. +type TestStatistic = Sample -> Sample -> Double + +-- | +-- Calculates the difference mean of two samples (/mean(s1 - s2)/). When the +-- two samples do not have an equal length, the trailing elements of the +-- longer vector are ignored. +differenceMean :: TestStatistic +differenceMean v1 v2 = + VG.sum (subVector v1 v2) / fromIntegral (VG.length v1) + +-- | Calculates the mean difference of two samples (/mean(s1) - mean(s2)/). +meanDifference :: TestStatistic +meanDifference s1 s2 = + mean s1 - mean s2 + +-- | Calculate the mean of a sample. +mean :: Sample -> Double +mean = do + t <- VG.sum + l <- VG.length + return $ t / fromIntegral l + +-- | Calculate the ratio of sample variances (/var(s1) : var(s2)/). +varianceRatio :: TestStatistic +varianceRatio v1 v2 = + variance v1 / variance v2 + +-- | Subtract two vectors. +subVector :: (VG.Vector v n, Num n) => v n -> v n -> v n +subVector = VG.zipWith (-) + +subIIW :: Int -> Int -> Word +subIIW a b = fromIntegral a - fromIntegral b +{-# INLINE subIIW #-} + +addIWI :: Int -> Word -> Int +addIWI a b = a + fromIntegral b +{-# INLINE addIWI #-} + +-- | Generate Int numbers within a range +randomIntR :: PureMT -> (Int, Int) -> (Int, PureMT) +randomIntR gen (a, b) + | n == 0 = randomInt gen + | otherwise = loop gen + where + (a', b') = if a < b then (a, b) else (b, a) + -- Number of different Ints that should be generated + n = 1 + subIIW b' a' + -- The total range of Word can hold x complete n ranges + x = maxBound `div` n + -- Pick from a range the is dividable by n without remainders + s = x * n + loop gen' + | r >= s = loop gen'' -- r is outside the range, discard it... + | otherwise = (addIWI a' (r `div` x), gen'') + where + (!r, !gen'') = randomWord gen' +{-# INLINE randomIntR #-} diff --git a/tests/tests.hs b/tests/tests.hs new file mode 100644 --- /dev/null +++ b/tests/tests.hs @@ -0,0 +1,53 @@ +module Main where + +import Control.Monad.Mersenne.Random (evalRandom) +import System.Random.Mersenne.Pure64 (pureMT) +import qualified Data.Vector.Unboxed as V +import Statistics.Test.ApproxRand +import Test.HUnit (assertEqual) +import Test.Framework +import Test.Framework.Providers.HUnit + +tests :: Test +tests = testGroup "Paired approximate randomization tests" $ + concat [statTests, randomizationTests] + +main :: IO () +main = defaultMain [ tests ] + +-- Statistics tests + +statTests :: [Test] +statTests = [meanDifferenceTest] + +meanDifferenceTest :: Test +meanDifferenceTest = + testEquality "mean difference robot competition" + (meanDifference cohenRobotsAlpha cohenRobotsBeta) 1.8 + +-- Approximate andomization tests + +randomizationTests :: [Test] +randomizationTests = [pairApproxExactTestScores] + +pairApproxExactTestScores :: Test +pairApproxExactTestScores = + testEquality "number of extreme values robot competition" + (length $ filter (>= 1.8) scores) 21 + where + scores = evalRandom + (approxRandPairScores meanDifference 1024 + cohenRobotsAlpha cohenRobotsBeta) $ + pureMT 42 + +-- Helper functions + +testEquality :: (Show a, Eq a) => String -> a -> a -> Test +testEquality msg a b = testCase msg $ assertEqual msg a b + +-- Example from Cohen, 1995 +cohenRobotsAlpha :: V.Vector Double +cohenRobotsAlpha = V.fromList [8,3,9,6,5,8,7,8,9,9] + +cohenRobotsBeta :: V.Vector Double +cohenRobotsBeta = V.fromList [7,0,9,4,5,9,8,3,4,5] diff --git a/utils/approx-rand-test-paired.hs b/utils/approx-rand-test-paired.hs new file mode 100644 --- /dev/null +++ b/utils/approx-rand-test-paired.hs @@ -0,0 +1,186 @@ +-- | +-- Copyright : (c) 2012 Daniël de Kok +-- License : Apache 2 +-- +-- Maintainer : Daniël de Kok +-- Stability : experimental +-- +-- Approximate randomization test (Noreen, 1989) + +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE DoAndIfThenElse #-} + +import Control.Exception.Base (Exception) +import Control.Monad (liftM, when) +import Control.Monad.Error (runErrorT) +import Control.Monad.Mersenne.Random (evalRandom) +import Data.Conduit (($$), ($=)) +import qualified Data.Conduit as C +import qualified Data.Conduit.Binary as CB +import qualified Data.Conduit.List as CL +import qualified Data.Conduit.Text as CT +import qualified Data.Text as T +import qualified Data.Text.Read as TR +import Data.Typeable (Typeable) +import qualified Data.Vector.Unboxed as V +import Data.Word (Word64) +import Statistics.Test.ApproxRand +import Statistics.Test.Types (TestType(..)) +import Statistics.Types (Sample) +import System.Console.GetOpt +import System.Environment (getArgs) +import System.Exit (exitFailure) +import System.Random.Mersenne.Pure64 (PureMT, newPureMT, pureMT) +import Text.Printf (printf) + +data ReadException = + DoubleConversionException String + deriving (Show, Typeable) + +instance Exception ReadException + +readFileCol :: String -> Int -> IO [Double] +readFileCol fn col = + liftM reverse $ C.runResourceT ( + CB.sourceFile fn $= + CB.lines $= + CT.decode CT.utf8 $= + CL.map (T.split (== ' ')) $= + CL.map (!! col) $= + toDouble $$ + CL.consume ) + +toDouble :: C.MonadThrow m => C.Conduit T.Text m Double +toDouble = CL.mapM $ \v -> + case TR.double v of + Left err -> C.monadThrow $ DoubleConversionException err + Right (d, _) -> return d + +main :: IO () +main = do + -- Read command-line options and arguments. + (opts, args) <- getOptions + + -- Read score files + let col = pred $ optColumn opts + v1 <- liftM V.fromList $ readFileCol (head args) col + v2 <- liftM V.fromList $ readFileCol (args !! 1) col + + let stat = optTestStatistic opts + + prng <- case optPRNGSeed opts of + Just seed -> return $ pureMT seed + Nothing -> newPureMT + + if optPrintScores opts then + printScores opts stat prng v1 v2 + else + applyTest opts stat prng v1 v2 + + +applyTest :: Options -> TestStatistic -> PureMT -> Sample -> + Sample -> IO () +applyTest opts stat prng v1 v2 = do + 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 testType = optTestType opts + + let pTest = optSigP opts + let pTail = case testType of + OneTailed -> pTest + TwoTailed -> pTest / 2 + + -- Test information + putStrLn $ "Test type: " ++ show testType + 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 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 + +printScores :: Options -> TestStatistic -> PureMT -> Sample -> + Sample -> IO () +printScores opts stat prng v1 v2 = do + let test = runErrorT $ approxRandPairScores stat (optIterations opts) v1 v2 + case evalRandom test prng of + Left err -> putStrLn err + Right scores -> mapM_ (putStrLn . printf "%f") scores + +data Options = Options { + optColumn :: Int, + optIterations :: Int, + optPRNGSeed :: Maybe Word64, + optPrintScores :: Bool, + optSigP :: Double, + optTestStatistic :: TestStatistic, + optTestType :: TestType +} + +defaultOptions :: Options +defaultOptions = Options { + optColumn = 1, + optIterations = 10000, + optPRNGSeed = Nothing, + optPrintScores = False, + optSigP = 0.01, + optTestStatistic = differenceMean, + optTestType = TwoTailed +} + +options :: [OptDescr (Options -> Options)] +options = + [ Option ['c'] ["column"] + (ReqArg (\arg opt -> opt { optColumn = read arg }) "NUMBER") + "column number (starting at 1)", + Option ['i'] ["iterations"] + (ReqArg (\arg opt -> opt { optIterations = read arg }) "NUMBER") + "number of iterations", + Option ['o'] ["one-tailed"] + (NoArg (\opt -> opt { optTestType = OneTailed })) + "perform a one-tailed test", + 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 ['s'] ["seed"] + (ReqArg (\arg opt -> opt { optPRNGSeed = Just $ read arg}) "NUMBER") + "pseudorandom number generator seed", + Option ['t'] ["test-statistic"] + (ReqArg (\arg opt -> opt { optTestStatistic = parseStatistic arg}) "NAME") + "test statistic (mean_diff, var_ratio)" + ] + +getOptions :: IO (Options, [String]) +getOptions = do + args <- getArgs + case getOpt Permute options args of + (actions, nonOpts, []) -> do + when (length nonOpts /= 2) usageExit + let opts = foldl (flip ($)) defaultOptions actions + return (opts, nonOpts) + (_, _, _) -> + usageExit + where + usageExit = do + putStrLn $ usageInfo header options + exitFailure + where + header = "Usage: approx-rand-test [OPTION...] scores scores2" + +parseStatistic :: String -> TestStatistic +parseStatistic "mean_diff" = differenceMean +parseStatistic "var_ratio" = varianceRatio +parseStatistic _ = error "Unknown test statistic" diff --git a/utils/approx-rand-test.hs b/utils/approx-rand-test.hs new file mode 100644 --- /dev/null +++ b/utils/approx-rand-test.hs @@ -0,0 +1,182 @@ +-- | +-- Copyright : (c) 2012 Daniël de Kok +-- License : Apache 2 +-- +-- Maintainer : Daniël de Kok +-- Stability : experimental +-- +-- Approximate randomization test (Noreen, 1989) + +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE DoAndIfThenElse #-} + +import Control.Exception.Base (Exception) +import Control.Monad (liftM, when) +import Control.Monad.Mersenne.Random (evalRandom) +import Data.Conduit (($$), ($=)) +import qualified Data.Conduit as C +import qualified Data.Conduit.Binary as CB +import qualified Data.Conduit.List as CL +import qualified Data.Conduit.Text as CT +import qualified Data.Text as T +import qualified Data.Text.Read as TR +import Data.Typeable (Typeable) +import qualified Data.Vector.Unboxed as V +import Data.Word (Word64) +import Statistics.Test.ApproxRand +import Statistics.Test.Types (TestType(..)) +import Statistics.Types (Sample) +import System.Console.GetOpt +import System.Environment (getArgs) +import System.Exit (exitFailure) +import System.Random.Mersenne.Pure64 (PureMT, newPureMT, pureMT) +import Text.Printf (printf) + +data ReadException = + DoubleConversionException String + deriving (Show, Typeable) + +instance Exception ReadException + +readFileCol :: String -> Int -> IO [Double] +readFileCol fn col = + liftM reverse $ C.runResourceT ( + CB.sourceFile fn $= + CB.lines $= + CT.decode CT.utf8 $= + CL.map (T.split (== ' ')) $= + CL.map (!! col) $= + toDouble $$ + CL.consume ) + +toDouble :: C.MonadThrow m => C.Conduit T.Text m Double +toDouble = CL.mapM $ \v -> + case TR.double v of + Left err -> C.monadThrow $ DoubleConversionException err + Right (d, _) -> return d + +main :: IO () +main = do + -- Read command-line options and arguments. + (opts, args) <- getOptions + + -- Read score files + let col = pred $ optColumn opts + v1 <- liftM V.fromList $ readFileCol (head args) col + v2 <- liftM V.fromList $ readFileCol (args !! 1) col + + let stat = optTestStatistic opts + + prng <- case optPRNGSeed opts of + Just seed -> return $ pureMT seed + Nothing -> newPureMT + + if optPrintScores opts then + printScores opts stat prng v1 v2 + else + applyTest opts stat prng v1 v2 + + +applyTest :: Options -> TestStatistic -> PureMT -> Sample -> + Sample -> IO () +applyTest opts stat prng v1 v2 = do + 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 testType = optTestType opts + + let pTest = optSigP opts + let pTail = case testType of + OneTailed -> pTest + TwoTailed -> pTest / 2 + + -- Test information + putStrLn $ "Test type: " ++ show testType + putStrLn $ printf "Test significance: %f" pTest + putStrLn $ printf "Tail significance: %f" pTail + + -- Approximate randomization testing. + let test = approxRandTest testType stat (optIterations opts) pTest v1 v2 + let result = evalRandom test prng + case result of + Significant p -> putStrLn $ printf "Significant: %f" p + NotSignificant p -> putStrLn $ printf "Not significant: %f" p + +printScores :: Options -> TestStatistic -> PureMT -> Sample -> + Sample -> IO () +printScores opts stat prng v1 v2 = + mapM_ (putStrLn . printf "%f") $ + evalRandom (approxRandScores stat (optIterations opts) v1 v2) prng + +data Options = Options { + optColumn :: Int, + optIterations :: Int, + optPRNGSeed :: Maybe Word64, + optPrintScores :: Bool, + optSigP :: Double, + optTestStatistic :: TestStatistic, + optTestType :: TestType +} + +defaultOptions :: Options +defaultOptions = Options { + optColumn = 1, + optIterations = 10000, + optPRNGSeed = Nothing, + optPrintScores = False, + optSigP = 0.01, + optTestStatistic = meanDifference, + optTestType = TwoTailed +} + +options :: [OptDescr (Options -> Options)] +options = + [ Option ['c'] ["column"] + (ReqArg (\arg opt -> opt { optColumn = read arg }) "NUMBER") + "column number (starting at 1)", + Option ['i'] ["iterations"] + (ReqArg (\arg opt -> opt { optIterations = read arg }) "NUMBER") + "number of iterations", + Option ['o'] ["one-tailed"] + (NoArg (\opt -> opt { optTestType = OneTailed })) + "perform a one-tailed test", + 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 ['s'] ["seed"] + (ReqArg (\arg opt -> opt { optPRNGSeed = Just $ read arg}) "NUMBER") + "pseudorandom number generator seed", + Option ['t'] ["test-statistic"] + (ReqArg (\arg opt -> opt { optTestStatistic = parseStatistic arg}) "NAME") + "test statistic (mean_diff, var_ratio)" + ] + +getOptions :: IO (Options, [String]) +getOptions = do + args <- getArgs + case getOpt Permute options args of + (actions, nonOpts, []) -> do + when (length nonOpts /= 2) usageExit + let opts = foldl (flip ($)) defaultOptions actions + return (opts, nonOpts) + (_, _, _) -> + usageExit + where + usageExit = do + putStrLn $ usageInfo header options + exitFailure + where + header = "Usage: approx-rand-test [OPTION...] scores scores2" + +parseStatistic :: String -> TestStatistic +parseStatistic "mean_diff" = meanDifference +parseStatistic "var_ratio" = varianceRatio +parseStatistic _ = error "Unknown test statistic"