Julia の HypothesisTests.PowerDivergenceTest() は問題がある。
対象とする分割表に，度数が 0 のセルがあると不適切な（役に立たない）結果を返す。
julia> x = [4 5 2 0; 0 7 6 1; 1 0 3 1]
3×4 Matrix{Int64}:
4 5 2 0
0 7 6 1
1 0 3 1
julia> using HypothesisTests
julia> result = PowerDivergenceTest(x, lambda=0.0);
julia> println(result.stat)
NaN
julia> println(pvalue(result))
NaN
そこで，以下のような関数を書く。
julia> function PDTsummary(result::PowerDivergenceTest; resid::Bool=false)
# 結果を簡潔に表示するための関数
function printmatrix(name, x)
println("\n$name")
Base.print_matrix(stdout, round.(x, digits=5), "", " ", "\n")
end
println("chisq. = $(round(result.stat, digits=5)), df = $(result.df), p.value = $(round(pvalue(result), digits=5))")
if resid
printmatrix("residuals", result.residuals)
printmatrix("standardized residuals", result.stdresiduals)
end
end
PDTsummary (generic function with 1 method)
julia> function g2stat(x, nrows, ncols, n, correct)
# セルに 0 があっても正しい答えを出す
# correct=true で，連続性の補正も行うことができる
ln(n) = n .== 0 ? 0 : n .* log.(n)
n = sum(x) # 全サンプルサイズ
n1 = sum(x, dims=2) # 行和
n2 = sum(x, dims=1) # 列和
G2 = 2*(sum(ln.(x)) - sum(ln.(n1)) - sum(ln.(n2)) + ln(n)) # G 統計量
correct && (G2 /= 1 + (n * sum(1 ./ n1) - 1) * (n * sum(1 ./ n2) - 1) / (6n * nrows * ncols)) # 連続性の補正
df = (nrows - 1) * (ncols - 1) # G の自由度
G2, df
end
g2stat (generic function with 1 method)
julia> function chisq_test(x::AbstractMatrix{T}; correct::Bool=false, resid::Bool=false, lambda::Float64=1.0) where T<:Integer
# 分割表を入力してχ二乗検定，G2 検定
nrows, ncols = size(x)
n = sum(x)
if lambda == 1.0 && correct && nrows == 2 && ncols == 2
a, c, b, d = x
stat = n*max(0, abs(a*d - b*c) - 0.5n)^2 / ((a+b)*(c+d)*(a+c)*(b+d))
result = PowerDivergenceTest(lambda, ones(4), stat, 1, x, n, ones(4), x, x, x)
elseif lambda == 1.0
result = PowerDivergenceTest(x; lambda)
elseif lambda == 0.0
G2, df = g2stat(x, nrows, ncols, n, correct)
vec = ones(nrows*ncols)
mat = reshape(vec, nrows, ncols)
result = PowerDivergenceTest(lambda, vec, G2, df, mat, n, vec, mat, mat, mat)
end
PDTsummary(result; resid)
end
chisq_test (generic function with 1 method)
julia> function G2_test(x::AbstractMatrix{T}; correct::Bool=false, resid::Bool=false) where T<:Integer
chisq_test(x; correct, resid, lambda=0.0)
end
G2_test (generic function with 1 method)
これで，
julia> G2_test(x)
chisq. = 15.36459, df = 6, p.value = 0.0176
という正しい結果を返す。
ついでに，連続性の補正もできるようになっているので，
julia> G2_test(x, correct=true)
chisq. = 13.77649, df = 6, p.value = 0.03224
のようになる。
