Repository: JuliaStats/HypothesisTests.jl
Pull request #147 (2019-01-26 07:47:55)

Add Barlett's test
Conflicts with #146 and more significantly with #145.

Imported from https://github.com/ararslan/MultivariateTests.jl, test values validated against SAS when originally written.


Patch:
diff --git a/docs/src/multivariate.md b/docs/src/multivariate.md
index ce7a0a2..cc8e670 100644
--- a/docs/src/multivariate.md
+++ b/docs/src/multivariate.md
@@ -7,3 +7,12 @@ OneSampleHotellingT2
 EqualCovHotellingT2
 UnequalCovHotellingT2
 ```
+
+## Equality of covariance matrices
+
+Bartlett's test for equality of two covariance matrices is provided.
+This is equivalent to Box's ``M``-test for two groups.
+
+```@docs
+BartlettsTest
+```
diff --git a/src/bartlett.jl b/src/bartlett.jl
new file mode 100644
index 0000000..7c8a332
--- /dev/null
+++ b/src/bartlett.jl
@@ -0,0 +1,60 @@
+# Tests for equality of covariance matrices
+
+export BartlettsTest
+
+abstract type CovarianceEqualityTest <: HypothesisTest end
+
+population_param_of_interest(::CovarianceEqualityTest) =
+    ("Equality of covariance matrices", NaN, NaN)
+
+## Utility functions
+
+# Finite population correction factor
+@inline _correction(p::Int, nxm1::Int, nym1::Int) =
+    1 - ((2p^2 + 3p - 1) * (inv(nxm1) + inv(nym1) - inv(nxm1 + nym1)) / 6(p + 1))
+
+## Bartlett's test
+
+struct BartlettsTest <: CovarianceEqualityTest
+    L′::Real
+    p::Int
+    nx::Int
+    ny::Int
+end
+
+"""
+    BartlettsTest(X::AbstractMatrix, Y::AbstractMatrix)
+
+Perform Bartlett's test of the hypothesis that the covariance matrices of `X` and `Y`
+are equal.
+
+!!! note
+    Bartlett's test is sensitive to departures from multivariate normality.
+"""
+function BartlettsTest(X::AbstractMatrix, Y::AbstractMatrix)
+    nx, p = size(X)
+    ny, q = size(Y)
+    p == q || throw(DimensionMismatch("Inconsistent number of variables"))
+    a = nx - 1
+    b = ny - 1
+    Sx = cov(X)
+    Sy = cov(Y)
+    L′ = -a * logdet(Sx) - b * logdet(Sy)
+    L′ += (a + b) * logdet(poolcov!(Sx, a, Sy, b))
+    L′ *= _correction(p, a, b)
+    return BartlettsTest(L′, p, nx, ny)
+end
+
+StatsBase.nobs(B::BartlettsTest) = (B.nx, B.ny)
+StatsBase.dof(B::BartlettsTest) = div(B.p * (B.p + 1), 2)
+
+testname(::BartlettsTest) = "Bartlett's Test for Equality of Covariance Matrices"
+default_tail(::BartlettsTest) = :right
+pvalue(B::BartlettsTest; tail=:right) = pvalue(Chisq(dof(B)), B.L′, tail=tail)
+
+function show_params(io::IO, B::BartlettsTest, indent="")
+    println(io, indent, "number of observations: ", nobs(B))
+    println(io, indent, "number of variables:    ", B.p)
+    println(io, indent, "χ² statistic:           ", B.L′)
+    println(io, indent, "degrees of freedom:     ", dof(B))
+end


Test patch:
diff --git a/src/HypothesisTests.jl b/src/HypothesisTests.jl
index 1e1d455..6dac75e 100644
--- a/src/HypothesisTests.jl
+++ b/src/HypothesisTests.jl
@@ -178,4 +178,5 @@ include("jarque_bera.jl")
 include("durbin_watson.jl")
 include("permutation.jl")
 include("hotelling.jl")
+include("bartlett.jl")
 end
diff --git a/test/bartlett.jl b/test/bartlett.jl
new file mode 100644
index 0000000..396bae0
--- /dev/null
+++ b/test/bartlett.jl
@@ -0,0 +1,18 @@
+using HypothesisTests
+using Test
+using StatsBase
+using DelimitedFiles
+
+@testset "Bartlett's test" begin
+    # Columns are type, length, left, right, bottom, top, diag
+    swiss = readdlm(joinpath(@__DIR__, "data", "swiss3.txt"))
+    genuine = convert(Matrix{Float64}, swiss[view(swiss, :, 1) .== "real", 2:end])
+    counterfeit = convert(Matrix{Float64}, swiss[view(swiss, :, 1) .== "fake", 2:end])
+
+    b = BartlettsTest(genuine, counterfeit)
+    @test nobs(b) == (100, 100)
+    @test dof(b) == 21
+    @test pvalue(b) ≈ 0.0 atol=1e-10
+    @test b.L′ ≈ 121.8991235 atol=1e-6
+    @test occursin("reject h_0", sprint(show, b))
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 2c6ecbd..dd216d1 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -32,3 +32,4 @@ include("t.jl")
 include("wilcoxon.jl")
 include("z.jl")
 include("hotelling.jl")
+include("bartlett.jl")