一元配置分散分析 なじみのない用語かも知れないが，「独立 k 標本の平均値の差の検定」である。 残念なことに，Julia では標準的に用意されておらず，GitHab には何例かあるものの，全ての人が簡単にアクセスできるものではないかもしれないので，あえてプログラムを書き，公開してみたい。 まず最初の関数は，測定データ x とそのデータがどの群に属するかを表す g で検定を行うという onewayAnova 関数である。 Julia> using Statistics, Distributions Julia> function onewayAnova(x, g; var_equal=false) Julia> member = Set(g) Julia> k = length(member) Julia> gi = [] Julia> ni = [] Julia> mi = [] Julia> vi = [] Julia> for i in member Julia> append!(gi, i) Julia> append!(ni, sum(g .== i)) Julia> append!(mi, mean(x[g .== i])) Julia> append!(vi, var(x[g .== i])) Julia> end Julia> method = "One-way analysis of means" Julia> if var_equal Julia> n = sum(ni) Julia> mx = mean(x) Julia> df1, df2 = k - 1, n - k Julia> F = ((sum(ni .* (mi .- mx).^2) / df1) / (sum((ni .- 1) .* vi) / df2)) Julia> else Julia> method *= " (not assuming equal variances)" Julia> wi = ni ./ vi Julia> sum_wi = sum(wi) Julia> tmp = sum((1 .- wi ./ sum_wi).^2 ./ (ni .- 1)) ./ (k^2 - 1) Julia> m = sum(wi .* mi) / sum_wi Julia> df1, df2 = k - 1, 1 / (3 * tmp) Julia> F = sum(wi .* (mi .- m).^2) / (df1 * (1 + 2 * (k - 2) * tmp)) Julia> end Julia> order = sortperm(gi) Julia> gi = gi[order] Julia> ni = ni[order] Julia> mi = mi[order] Julia> vi = vi[order] Julia> output1(gi, ni, mi, vi) Julia> output2(method, F, df1, df2) Julia> end onewayAnova (generic function with 2 methods) もう一つの関数は，なんと，同じ名前 onewayAnova であるが，引数（の型と個数）が違う。 Julia は，同じ名前の関数があると，呼ばれたのはどの関数だろう？と判断して，適切な関数を呼んでくれる。 なんと，賢い！！ Julia> function onewayAnova(ni, mi, ui, gi; var_equal=false) Julia> x = [] Julia> g = [] Julia> for i in 1:length(ni) Julia> obj = Normal(mi[i], sqrt(ui[i])) Julia> temp = rand(obj, ni[i]) Julia> temp = (temp .- mean(temp)) ./ std(temp) .* sqrt(ui[i]) .+ mi[i] Julia> append!(x, temp) Julia> append!(g, repeat([gi[i]], ni[i])) Julia> end Julia> ve = var_equal Julia> onewayAnova(x, g, var_equal=ve) Julia> end onewayAnova (generic function with 2 methods) 次の 2 つの関数は，検定の結果を出力するための関数である。 Julia> using Printf Julia> function output1(gi, ni, mi,vi) Julia> println("Sumary of the data.") Julia> @printf("%10s: %5s %15s %15s %15s\n", "group", "n", "mean", "variance", "s.d.") Julia> for i = 1:length(gi) Julia> @printf("%10s: %5d %15.7f %15.7f %15.7f\n", gi[i], ni[i], mi[i], vi[i], sqrt(vi[i])) Julia> end Julia> end Julia> function output2(method, F, df1, df2) Julia> println(method) Julia> obj = FDist(df1, df2) Julia> if df2 == floor(df2) Julia> @printf("F = %.7g, d.f.1 = %d, d.f.2 = %d, p value = %.7g\n", F, df1, df2, 1 - cdf(obj, F)) Julia> else Julia> @printf("F = %.7g, d.f.1 = %d, d.f.2 = %.7g, p value = %.7g\n", F, df1, df2, 1 - cdf(obj, F)) Julia> end Julia> end output2 (generic function with 1 method) 使用例 最初の使用例は，データベクトルと，群データベクトルを用いて呼び出すものである。 まず，R ではデフォルトではないが，群の分散が等しいと仮定する場合である。 Julia> # > x = c(2, 1, 2, 3, 4, 3, 5, 4, 6, 7, 2, 4) Julia> # > g = c("a", "a", "a", "a", "b", "c", "c", "b", "b", "c", "c", "c") Julia> # > oneway.test(x ~ g, var.equal=TRUE) Julia> # Julia> # One-way analysis of means Julia> # Julia> # data: x and g Julia> # F = 3.5715, num df = 2, denom df = 9, p-value = 0.07214 Julia> ##### Julia> x = [2, 1, 2, 3, 4, 3, 5, 4, 6, 7, 2, 4] Julia> g = ["a", "a", "a", "a", "b", "c", "c", "b", "b", "c", "c", "c"] Julia> onewayAnova(x, g, var_equal=true) Sumary of the data. group: n mean variance s.d. a: 4 2.0000000 0.6666667 0.8164966 b: 3 4.6666667 1.3333333 1.1547005 c: 5 4.2000000 3.7000000 1.9235384 One-way analysis of means F = 3.57149, d.f.1 = 2, d.f.2 = 9, p value = 0.0721381 次の使用例は，同じデータに対して，群の分散が等しいことを仮定しない場合（独立 2 標本の t 検定の場合の，Welch の方法の拡張）である。 これは R ではデフォルトになっている。 Julia> # > oneway.test(x ~ g) Julia> # Julia> # One-way analysis of means (not Julia> # assuming equal variances) Julia> # Julia> # data: x and g Julia> # F = 6.2568, num df = 2.0000, denom df = 5.0833, p-value = 0.04259 Julia> ##### Julia> onewayAnova(x, g) Sumary of the data. group: n mean variance s.d. a: 4 2.0000000 0.6666667 0.8164966 b: 3 4.6666667 1.3333333 1.1547005 c: 5 4.2000000 3.7000000 1.9235384 One-way analysis of means (not assuming equal variances) F = 6.256838, d.f.1 = 2, d.f.2 = 5.083311, p value = 0.04258988 二番目の例は，群ごとに，「サンプルサイズ」，「平均値」，「不偏分散」，「群の識別子」を与えて検定を行う方法である。 Julia> onewayAnova([4, 3, 5], [2, 4.6666667, 4.2], [0.6666667, 1.3333333, 3.7], ["a", "b", "c"], var_equal=true) Sumary of the data. group: n mean variance s.d. a: 4 2.0000000 0.6666667 0.8164966 b: 3 4.6666667 1.3333333 1.1547005 c: 5 4.2000000 3.7000000 1.9235384 One-way analysis of means F = 3.57149, d.f.1 = 2, d.f.2 = 9, p value = 0.07213809 同じデータで，等分散を仮定しない場合である。 Julia> onewayAnova([4, 3, 5], [2, 4.6666667, 4.2], [0.6666667, 1.3333333, 3.7], ["a", "b", "c"]) Sumary of the data. group: n mean variance s.d. a: 4 2.0000000 0.6666667 0.8164966 b: 3 4.6666667 1.3333333 1.1547005 c: 5 4.2000000 3.7000000 1.9235384 One-way analysis of means (not assuming equal variances) F = 6.256838, d.f.1 = 2, d.f.2 = 5.083311, p value = 0.04258988
1. Julia の F 分布関数
Rmath パッケージを使用すれば，R の pf(), qf() が使える。第 4 引数は lower.tail に対応し，true （デフォルト）のとき lower.tail=TRUE，false のとき lower.tail=FALSE を意味する。
using Distributions, Rmath
1.1. オブジェクト生成 FDist( )
obj = FDist(3, 10)
FDist{Float64}(ν1=3.0, ν2=10.0)
1.2. 上側確率 ccdf( ), pf( )
ccdf(obj, 1.5)
0.2737765559785967
pf(1.5, 3, 10, false)
0.2737765559785967
1.3. 上側確率に対するパーセント点 cquantile( ), qf( )
cquantile(obj, 0.05)
3.7082648190468457
qf(0.05, 3, 10, false)
3.7082648190468457
1. Julia の t 分布関数
Rmath パッケージを使用すれば，R の pt(), qt() が使える。第 3 引数は lower.tail に対応し，true （デフォルト）のとき lower.tail=TRUE，false のとき lower.tail=FALSE を意味する。
using Distributions, Rmath
1.1. オブジェクト生成 TDist(ν)
obj10 = TDist(10)
TDist{Float64}(ν=10.0)
1.2. 確率
1.2.1. 上側確率 ccdf( ), pt( )
ccdf(obj10, 1.96)
0.03921812012384987
pt(1.96, 10, false)
0.03921812012384987
1.2.2. Julia 下側確率 cdf( ), pt( )
cdf(obj10, -1.96)
0.03921812012384987
pt(-1.96, 10)
0.03921812012384987
1.3. パーセント点
1.3.1. 上側確率に対するパーセント点 cquantile( ), qt( )
cquantile(obj10, 0.025)
2.2281388519862744
qt(0.025, 10, false)
2.2281388519862744
1.3.2. 下側確率に対するパーセント点 quantile( ), qt( )
quantile(obj10, 0.025)
-2.2281388519862744
qt(0.025, 10)
-2.2281388519862744
1. Julia のカイ二乗分布関数
Rmath パッケージを使用すれば，R の pchisqt(), qchisq() が使える。第 2 引数は lower.tail に対応し，true （デフォルト）のとき lower.tail=TRUE，false のとき lower.tail=FALSE を意味する。
using Distributions, Rmath
1.1. オブジェクト生成
obj1 = Chisq(1)
Chisq{Float64}(ν=1.0)
1.2. 上側確率
ccdf(obj1, 3.841459)
0.04999999465319573
pchisq(3.841459, 1, false)
0.04999999465319573
1.3. 上側確率に対するパーセント点
cquantile(obj1, 0.05)
3.841458820694126
qchisq(0.05, 1, false)
3.841458820694126
PVアクセスランキング にほんブログ村
