packages feed

approx-rand-test (empty) → 0.0.3

raw patch · 7 files changed

+1040/−0 lines, 7 filesdep +HUnitdep +approx-rand-testdep +basesetup-changed

Dependencies added: HUnit, approx-rand-test, base, conduit, mersenne-random-pure64, monad-mersenne-random, mtl, resourcet, statistics, test-framework, test-framework-hunit, text, transformers, vector

Files

+ LICENSE view
@@ -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.
+ Setup.lhs view
@@ -0,0 +1,4 @@+#!/usr/bin/env runhaskell++> import Distribution.Simple+> main = defaultMain
+ approx-rand-test.cabal view
@@ -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 <me@danieldk.eu>+Maintainer:    Daniël de Kok <me@danieldk.eu>+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.*
+ src/Statistics/Test/ApproxRand.hs view
@@ -0,0 +1,342 @@+-- |+-- Copyright  : (c) 2012 Daniël de Kok+-- License    : BSD3+--+-- Maintainer : Daniël de Kok <me@danieldk.eu>+-- 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 #-}
+ tests/tests.hs view
@@ -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]
+ utils/approx-rand-test-paired.hs view
@@ -0,0 +1,186 @@+-- |+-- Copyright  : (c) 2012 Daniël de Kok+-- License    : Apache 2+--+-- Maintainer : Daniël de Kok <me@danieldk.eu>+-- 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"
+ utils/approx-rand-test.hs view
@@ -0,0 +1,182 @@+-- |+-- Copyright  : (c) 2012 Daniël de Kok+-- License    : Apache 2+--+-- Maintainer : Daniël de Kok <me@danieldk.eu>+-- 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"