Merge pull request #76 from bjarthur/permutationtest

added permutation test

docs/src/nonparametric.md

@@ -59,3 +59,10 @@ SignedRankTest
 ExactSignedRankTest
 ApproximateSignedRankTest
 ```
+
+## Permutation test
+
+```@docs
+ExactPermutationTest
+ApproximatePermutationTest
+```

src/HypothesisTests.jl

@@ -27,7 +27,7 @@ __precompile__()
 module HypothesisTests
 
 using Distributions, Roots, StatsBase, Compat
-using Combinatorics: combinations
+using Combinatorics: combinations, permutations
 using Rmath: pwilcox, psignrank
 
 import StatsBase.confint
@@ -175,4 +175,5 @@ include("breusch_godfrey.jl")
 include("augmented_dickey_fuller.jl")
 include("jarque_bera.jl")
 include("durbin_watson.jl")
+include("permutation.jl")
 end

src/permutation.jl

@@ -0,0 +1,59 @@
+export ExactPermutationTest, ApproximatePermutationTest
+
+function ptstats(x,y,f)
+    xy = vcat(x,y)
+    rx = 1:length(x)
+    ry = (length(xy)-length(y)+1):length(xy)
+    (xy, rx, ry)
+end
+
+struct PermutationTest{T<:Real} <: HypothesisTest
+    observation::T
+    samples::Vector{T}
+end
+
+"""
+    ExactPermutationTest(x::Vector, y::Vector, f::Function)
+
+Perform a permutation test (a.k.a. randomization test) of the null hypothesis
+that `f(x)` is equal to `f(y)`.  All possible permutations are sampled.
+"""
+function ExactPermutationTest(x::AbstractVector{R}, y::AbstractVector{S},
+                              f::Function) where {R<:Real,S<:Real}
+    xy, rx, ry = ptstats(x,y,f)
+    P = permutations(xy)
+    samples = [f(view(p,rx)) - f(view(p,ry)) for p in P]
+    PermutationTest(f(x) - f(y), samples)
+end
+
+"""
+    ApproximatePermutationTest(x::Vector, y::Vector, f::Function, n::Int)
+
+Perform a permutation test (a.k.a. randomization test) of the null hypothesis
+that `f(x)` is equal to `f(y)`.  `n` of the `factorial(length(x)+length(y))`
+permutations are sampled at random.
+"""
+function ApproximatePermutationTest(x::AbstractVector{R}, y::AbstractVector{S},
+                                    f::Function, n::Int) where {R<:Real,S<:Real}
+    xy, rx, ry = ptstats(x,y,f)
+    samples = [(shuffle!(xy); f(view(xy,rx)) - f(view(xy,ry))) for i = 1:n]
+    PermutationTest(f(x) - f(y), samples)
+end
+
+function pvalue(apt::PermutationTest; tail=:both)
+    if tail == :both
+        count = sum(abs(apt.observation) >= abs(x) for x in apt.samples)
+    elseif tail == :left
+        count = sum(apt.observation >= x for x in apt.samples)
+    elseif tail == :right
+        count = sum(apt.observation <= x for x in apt.samples)
+    end
+    return count / length(apt.samples)
+end
+
+testname(::PermutationTest) = "Permutation Test"
+
+function show_params(io::IO, apt::PermutationTest, ident)
+    println(io, ident, "observation: ", apt.observation)
+    println(io, ident, "samples: ", apt.samples)
+end

test/permutation.jl

@@ -0,0 +1,12 @@
+using HypothesisTests, Base.Test
+
+@test isapprox(pvalue(ExactPermutationTest([1,2,3],[4,5,6],mean), tail=:both), 1.0)
+@test isapprox(pvalue(ExactPermutationTest([4,5,6],[1,2,3],mean), tail=:both), 1.0)
+@test isapprox(pvalue(ExactPermutationTest([1,2,3],[4,5,6],mean), tail=:left), 0.05)
+@test isapprox(pvalue(ExactPermutationTest([4,5,6],[1,2,3],mean), tail=:right), 0.05)
+
+srand(12345)
+@test isapprox(pvalue(ApproximatePermutationTest([1,2,3],[4,5,6],mean,100), tail=:both), 1.0)
+@test isapprox(pvalue(ApproximatePermutationTest([4,5,6],[1,2,3],mean,100), tail=:both), 1.0)
+@test isapprox(pvalue(ApproximatePermutationTest([1,2,3],[4,5,6],mean,100), tail=:left), 0.08)
+@test isapprox(pvalue(ApproximatePermutationTest([4,5,6],[1,2,3],mean,100), tail=:right), 0.07)

test/runtests.jl

@@ -23,6 +23,7 @@ include("jarque_bera.jl")
 include("kolmogorov_smirnov.jl")
 include("kruskal_wallis.jl")
 include("mann_whitney.jl")
+include("permutation.jl")
 include("power_divergence.jl")
 include("t.jl")
 include("wilcoxon.jl")