using HypothesisTests
★ 独立性の検定
ChisqTest() は PowerDIvergenceTest() で，lambda=1.0 としたものである。
実際，HypothesisTests.jl で
function ChisqTest(x::AbstractMatrix{T}) where T<:Integer
PowerDivergenceTest(x, lambda=1.0)
end
と定義されている。
最も単純な 2×2 分割表での独立性の検定を行ってみる。
x = [11 6; 7 14]
#=
2×2 Matrix{Int64}:
11 6
7 14
=#
この分割表は，λの値により p value = 0.05 より大きかったり小さかったりする，微妙な領域にある例である（最後の方の図を参照）。
検定統計量は以下の式で計算される。
☆ λ = 1 の場合
PowerDIvergenceTest() と ChisqTest() の結果が同じであることを確認する。
統計量などは obj.stat，p 値は pvalue(obj) のように取り出せる。
println(obj) などとすると，全ての統計量が表示される。
obj1 = PowerDivergenceTest(x); # lambda=1.0 はデフォルト
println("statistic: $(obj1.stat), p value = $(pvalue(obj1))")
# statistic: 3.708932461873637, p value = 0.05412199588766715
obj2 = ChisqTest(x);
println("statistic: $(obj2.stat), p value = $(pvalue(obj2))")
# statistic: 3.708932461873637, p value = 0.05412199588766715
☆ λ の設定により，様々な検定統計量を求めることができる
☆ λ = 1 のとき，ピアソンのχ2統計量になる
obj3 = PowerDivergenceTest(x, lambda=1.0);
println("statistic: $(obj3.stat), p value = $(pvalue(obj3))")
# statistic: 3.708932461873637, p value = 0.05412199588766715
https://blog.goo.ne.jp/r-de-r/e/03fdab8cb34bab738b24abdc05117f35
☆ λ → 0 のとき，尤度比検定統計量（G2）に収束する
obj4 = PowerDivergenceTest(x, lambda=0.0);
println("statistic: $(obj4.stat), p value = $(pvalue(obj4))")
# statistic: 3.7658347787910724, p value = 0.052309741049692396
# 尤度比検定 G2
using Rmath
function g2(mat)
ln(n) = sum([x * log(x) for x in n if x > 0])
n = sum(mat)
rowsums = sum(mat, dims=2)
colsums = sum(mat, dims=1)
G2 = 2 * (ln(mat) - ln(rowsums) - ln(colsums) + ln(n))
nrow, ncol = size(mat)
df = (nrow - 1) * (ncol - 1)
P = pchisq(G2, df, false)
println("G2 = $G2 df = $df p value = $P")
end
g2(x) # G2 = 3.7658347787910884 df = 1 p value = 0.052309741049691945
https://blog.goo.ne.jp/r-de-r/e/4c8a247c619439dca239b1bc790bbf98
☆ λ → −1 のとき，最小判別情報統計量に収束する(Gokhale and Kullback, 1978)
obj5 = PowerDivergenceTest(x, lambda=-0.9999999999);
println("statistic: $(obj5.stat), p value = $(pvalue(obj5))")
# statistic: 3.89311661820617, p value = 0.04848437211114625
The minimum discrimination information approach in analyzing categorical data
D.V Gokhale & S. Kullback
Pages 987-1005 | Received 01 Nov 1977, Published online: 27 Jun 2007
☆ λ = −2 のとき，ネイマンの修正χ2に等しい(Neyman, 1949)
obj6 = PowerDivergenceTest(x, lambda=-2.0);
println("statistic: $(obj6.stat), p value = $(pvalue(obj6))")
# statistic: 4.0990514563921785, p value = 0.04290727930312255
Neyman, J. (1949) Contribution to the Theory of Chi-Square Test. Proceedings of the First Berkeley Symposium on Mathematical Statistics and Probability, 239-273.
☆ λ = −1/2 のとき，フリーマン・チューキー統計量である(Freeman and Tukey, 1950).
obj7 = PowerDivergenceTest(x, lambda=-0.5);
println("statistic: $(obj7.stat), p value = $(pvalue(obj7))")
# statistic: 3.820238526805399, p value = 0.05063702908109902
Transformations Related to the Angular and the Square Root
Murray F. Freeman, John W. Tukey
Ann. Math. Statist. 21(4): 607-611 (December, 1950)
☆ λ と検定統計量 stat と p value のグラフ化
x = [11 6; 7 14]
stat = []
pvalues = []
lambda = 1:-0.01:-2
for λ in lambda
obj = PowerDivergenceTest(x, lambda=λ)
append!(stat, obj.stat)
append!(pvalues, pvalue(obj))
end
using Plots
pyplot()
p1 = plot(lambda, stat, label="", tick_direction=:out,
xlabel="λ", ylabel="statistic");
p2 = plot(lambda, pvalues, label="", tick_direction=:out,
xlabel="λ", ylabel="p value");
plot(p1, p2, size=(400, 200))
savefig("fig.png")
有意になりやすいからといって lambda=-2 を指定するというのは邪道。
