Julia に翻訳--214 二群の平均値の差の検定

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

二群の平均値の差の検定
http://aoki2.si.gunma-u.ac.jp/R/my-t-test.html

ファイル名: myttest.jl
関数名:   myttest(x, y; varequal = false)
         myttest(nx, mx, ux, ny, my, uy; varequal = false)

翻訳するときに書いたメモ

==========#

using Statistics, Rmath, Printf

function myttest(x, y; varequal = false)
    myttest(length(x), mean(x), var(x), length(y), mean(y), var(y); varequal)
end

function myttest(nx, mx, ux, ny, my, uy; varequal = false)
    if varequal
        df = nx + ny - 2
        v = ((nx - 1) * ux + (ny - 1) * uy) / df
        tstat = abs(mx - my) / sqrt(v * (1 / nx + 1 / ny))
    else
        tstat = abs(mx - my) / sqrt(ux / nx + uy / ny)
        df = (ux / nx + uy / ny) ^ 2 / ((ux / nx) ^ 2 / (nx - 1) + (uy / ny) ^ 2 / (ny - 1))
    end
    pvalue = 2pt(tstat, df, false)
    @printf("t = %.5f,  df = %.5f,  p value = %.5f\n", tstat, df, pvalue)
end

x = [44, 50, 52, 53, 49, 53, 47, 44, 38, 62];
y = [60, 55, 49, 58, 72, 69, 61, 63, 61, 55, 59, 63];
myttest(x, y, varequal=true) # t = 4.13039,  df = 20.00000,  p value = 0.00052
myttest(x, y)                # t = 4.10725,  df = 18.82177,  p value = 0.00061

myttest(36, 82.6, 15.3, 43, 84.5, 16.2, varequal=true) # t = 2.11652,  df = 77.00000,  p value = 0.03753
myttest(36, 82.6, 15.3, 43, 84.5, 16.2)                # t = 2.12195,  df = 75.26729,  p value = 0.03713

Julia に翻訳--213 母平均の検定と推定

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

母平均の検定と推定
http://aoki2.si.gunma-u.ac.jp/lecture/Average/Mean1-r.html

ファイル名: onesamplettest.jl
関数名:   onesamplettest(x, mu0; conflevel=0.95)
         onesamplettest(n, xbar, u, mu0; conflevel=0.95)

翻訳するときに書いたメモ

==========#

using Statistics, StatsBase, Rmath, Printf

function onesamplettest(x, mu0; conflevel=0.95)
    n = length(x)
    xbar = mean(x)
    u = var(x)
    se = sqrt(u/n)
    t0 = abs(xbar - mu0) / se
    df = n - 1
    pvalue = 2pt(t0, df, false)
    @printf("t = %.5f,  df = %d,  p value = %.5f\n", t0, df, pvalue)
    t1 = qt((1 + conflevel)/2, df)
    lcl, ucl = xbar .+ [-1, 1]t1 .* se
    @printf("%s%% confidence interval = [%.5f, %.5f]\n", 100conflevel, lcl, ucl)
end

function onesamplettest(n, xbar, u, mu0; conflevel=0.95)
    x = randn(n);
    x = (x .- mean(x)) ./ std(x) .* sqrt(u) .+ xbar;
    onesamplettest(x, mu0; conflevel)
end

onesamplettest(31, 157.8, 24.6, 156.2)
# t = 1.79611,  df = 30,  p value = 0.08255
# 95.0% confidence interval = [155.98072, 159.61928]

x = [9.7, 10.3, 9.6, 7.7, 10.2, 10.6, 10.4, 11.4, 7.8, 8.6];
onesamplettest(x, 10.0)
# t = 0.95248,  df = 9,  p value = 0.36573
# 95.0% confidence interval = [8.75125, 10.50875]

onesamplettest(length(x), mean(x), var(x), 10)
# t = 0.95248,  df = 9,  p value = 0.36573
# 95.0% confidence interval = [8.75125, 10.50875]

Julia に翻訳--212 二群の比率の差の検定

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

二群の比率の差の検定
http://aoki2.si.gunma-u.ac.jp/lecture/Hiritu/diff-p-test-r.html

ファイル名: proptest.jl
関数名:    proptest(x::Vector{Int64}, n::Vector{Int64}; correction=true)
          proptest(x::Array{Int64,2}; correction=true)

翻訳するときに書いたメモ

K 群の比率の差の検定は，2xK 分割表での独立性の検定と同じである。ということで，chisqtest2 関数を下請けにしている。
下のプログラムは，2群の比率の差の検定だけでなく，K群の比率の差の検定にも対応している。

==========#

using Rmath, Printf

function proptest(x::Vector{Int64}, n::Vector{Int64}; correction=true)
    chisqtest2(hcat(x, n .- x); correction)
end

function proptest(x::Array{Int64,2}; correction=true)
    chisqtest2(x; correction)
end

function chisqtest2(x; correction=true)
    n, m = size(x)
    rowsum = sum(x, dims=2)
    colsum = sum(x, dims=1)
    expectation =  (rowsum * colsum) ./ sum(x)
    if sum(expectation .< 1) >= 1 ||sum(expectation .< 5) >= 0.2n*m
        println("Warning message: Chi-squared approximation may be incorrect")
    end
    difference = abs.(x - expectation)
    yates = 0
    correction && n == 2 && m == 2 && (yates = min(0.5, minimum(difference)))
    chisq = sum(((difference .- yates) .^ 2) ./ expectation)
    df =(n - 1) * (m - 1)
    pvalue = pchisq(chisq, df, false)
    @printf("chisq = %.5f,  df = %d,  p value = %.5f", chisq, df, pvalue)
end

smokers  = [83, 90, 129, 70]
patients = [86, 93, 136, 82];
proptest(smokers, patients) # chisq = 12.60041,  df = 3,  p value = 0.00559

data = [83 90 129 70
         3  3   7 12];
proptest(data)              # chisq = 12.60041,  df = 3,  p value = 0.00559

proptest([145, 157], [300, 250], correction=false) # chisq = 11.52663,  df = 1,  p value = 0.00069
proptest([145, 157], [300, 250])                   # chisq = 10.94973,  df = 1,  p value = 0.00094

Julia に翻訳--211 分布の差の検定，独立性の検定，K×M 分割表，2×2 分割表

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

分布の差の検定（独立性の検定）
http://aoki2.si.gunma-u.ac.jp/lecture/Cross/differenceofdist-r.html

独立性の検定（K×M 分割表）
http://aoki2.si.gunma-u.ac.jp/lecture/Cross/cross-r.html

独立性の検定（2×2 分割表）
http://aoki2.si.gunma-u.ac.jp/lecture/Cross/cross-r2.html

ファイル名: chisqtest2.jl　　関数名: chisqtest2

翻訳するときに書いたメモ

==========#

using Rmath, Printf

function chisqtest2(x; correction=true)
    n, m = size(x)
    rowsum = sum(x, dims=2)
    colsum = sum(x, dims=1)
    expectation =  (rowsum * colsum) ./ sum(x)
    if sum(expectation .< 1) >= 1 ||sum(expectation .< 5) >= 0.2n*m
        println("Warning message: Chi-squared approximation may be incorrect")
    end
    difference = abs.(x - expectation)
    yates = 0
    correction && n == 2 && m == 2 && (yates = min(0.5, minimum(difference)))
    chisq = sum(((difference .- yates) .^ 2) ./ expectation)
    df =(n - 1) * (m - 1)
    pvalue = pchisq(chisq, df, false)
    @printf("chisq = %.5f,  df = %d,  p value = %.5f", chisq, df, pvalue)
end

x = [20 15 16 4
     15  7  9 4];

chisqtest2(x)
# Warning message: Chi-squared approximation may be incorrect
# chisq = 1.19810,  df = 3,  p value = 0.75346

x = [4 2
     1 6];

chisqtest2(x, correction=false)
# Warning message: Chi-squared approximation may be incorrect
# chisq = 3.74524,  df = 1,  p value = 0.05296

chisqtest2(x)
# Warning message: Chi-squared approximation may be incorrect
# chisq = 1.85908,  df = 1,  p value = 0.17273

Julia に翻訳--210 適合度の検定，名義尺度の場合，χ2 分布による検定，一様性の検定，理論比が与えられる場合

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

適合度の検定－－名義尺度の場合；χ2 分布による検定；一様性の検定
http://aoki2.si.gunma-u.ac.jp/lecture/GoodnessOfFitness/nominalscale-r.html

適合度の検定－－名義尺度の場合；χ2 分布による検定；理論比が与えられる場合
http://aoki2.si.gunma-u.ac.jp/lecture/GoodnessOfFitness/nominalscale-r2.html

ファイル名: chisqtest.jl　　関数名: chisqtest

翻訳するときに書いたメモ

==========#

using Rmath, Printf

function chisqtest(x; p = fill(1 / length(x), length(x)))
    1 - sum(p) > eps() && error("sum(p) isn't equal to 1")
    expectation =  sum(x) .* p
    chisq = sum(((x .- expectation) .^ 2) ./ expectation)
    df = length(x) - 1
    pvalue = pchisq(chisq, df, false)
    @printf("chisq = %.5f,  df = %d,  p value = %.5f", chisq, df, pvalue)
end

x = [10, 12, 9, 4, 13, 8]
chisqtest(x)                         # chisq = 5.50000,  df = 5,  p value = 0.35795
chisqtest(x, p = fill(1/6, 6))       # chisq = 5.50000,  df = 5,  p value = 0.35795

y = [29, 12, 8, 2]
chisqtest(y, p = repeat([0.25], 4))  # chisq = 31.58824,  df = 3,  p value = 0.00000
chisqtest(y, p = [9, 3, 3, 1] ./ 16) # chisq = 1.32244,  df = 3,  p value = 0.72381