Julia 脳

Julia を始めたのが 2020 年の 12 月半ば。

もう，丸 7 ヶ月を過ぎようとしている。

R や Python のプログラムを書こうとしたら，「あ？こんなとき R （或いは Python）で，どう書くんだっけ？」と思うのはまだしも，Julia のままで無意識に書いて，シンタックスエラーが出まくりということになっている。

まあ，言語を切り替えるときにはよくあることではあるが。

で，ようやく？Julia で何の不都合もないところまで来た。思えば遠くへきたもんだ。（by 武田鉄矢）

Julia で正規性の検定，Anderson-Darling normality test, R nortest::ad.test

# https://www.spcforexcel.com/knowledge/basic-statistics/anderson-darling-test-for-normality
# https://www.spcforexcel.com/knowledge/basic-statistics/normal-probability-plots
# https://www.rdocumentation.org/packages/nortest/versions/1.0-4/topics/ad.test

R の nortest::ad.test では，修正前の検定統計量が返されるが，p 値は修正後の検定統計量に基づいたものである。

using Statistics, Rmath, Plots

function AndersonDarling(x; plot=true, color=:blue, alpha=0.4,
            xlabel="Sorted,Data", ylabel="Z value",
            title="Normal Probability Plot")
    n = length(sort!(x))
    Fx = pnorm.(x, mean(x), std(x));
    s = sum((2i - 1)*(log(Fx[i]) + log(1 - Fx[n-i+1])) for i = 1:n)
    AD = -n - s/n
    ADast = ((2.25/n + 0.75)/n + 1)AD
    if ADast >= 10
        pvalue = 3.7e-24
    elseif ADast >= 0.6
        pvalue = exp((0.0186ADast - 5.709)ADast + 1.2937)
    elseif ADast >= 0.34
        pvalue = exp((-1.38ADast - 4.279)ADast + 0.9177)
    elseif ADast >= 0.2
        pvalue = 1 - exp((-59.938ADast + 42.796)ADast - 8.318)
    else
        pvalue = 1 - exp((-223.73ADast + 101.14)ADast - 13.436)
    end
    if plot == true
        pyplot()
        z = qnorm.(((1:n) .- 0.5) / n)
        plt = scatter(x, z, grid=false, tick_direction=:out,
            color=color, markerstrokewidth=0, alpha=alpha,
            smooth=true,
            xlabel=xlabel, ylabel=ylabel, title=title, label="")
        display(plt)
    end
    println("AD = $AD,  AD* = $ADast,  p value = $pvalue")
    (AD=AD, ADast=ADast, pvalue=pvalue)
end

使用例

ForearmLengths = [
    17.3,20.9,18.7,17.9,18.3,19.0,18.1,18.8,19.1,17.9,
    18.2,19.4,19.4,17.3,18.3,19.0,20.5,18.5,19.4,19.6,
    19.0,20.4,18.6,18.3,19.6,20.4,16.1,19.6,19.3,21.0,
    18.3,18.7,18.5,17.2,18.0,19.9,18.8,20.0,17.5,17.9,
    18.7,17.3,17.8,19.6,18.1,20.9,18.1,19.8,17.6,19.5,
    17.7,19.9,16.6,20.0,17.1,19.1,19.6,19.4,19.9,18.9,
    19.7,18.4,19.3,16.9,18.5,18.1,19.5,20.1,19.5,19.2,
    18.4,16.8,20.5,20.4,20.5,17.5,17.1,20.0,19.1,18.3,
    18.9,18.9,20.8,18.5,19.4,19.0,19.7,17.7,18.3,21.4,
    20.5,19.7,19.9,19.8,19.0,17.3,19.2,18.8,19.1,18.6,
    18.3,20.6,16.4,17.5,19.5,18.4,20.1,18.5,18.5,17.4,
    18.6,18.8,19.0,19.3,18.5,19.8,17.1,20.6,19.1,18.4,
    20.2,18.6,19.2,17.4,18.3,18.5,18.0,17.1,16.3,20.7,
    18.5,18.7,16.3,18.2,19.3,18.0,20.3,17.2,18.8,17.7];

AndersonDarling(ForearmLengths, xlabel="Forearm Lengths");
# AD = 0.23699900708288624,  AD* = 0.2382958511395005,  p value = 0.7820446056638978
savefig("fig1.png")

R での結果

using RCall
R"""
library(nortest)
options(digits=15)
a = ad.test($ForearmLengths)
cat(a$statistic, a$p.value, "\n")
"""
# 0.236999007082858 0.782044605663986 

BirthWeights = [
    3837,3480,3334,3116,3554,3428,3838,3783,3625,
    3345,2208,3034,1745,2184,2846,3300,3166,2383,
    3520,3428,3380,4162,3294,3630,2576,3406,3208,
    3402,3521,3500,3746,3736,3523,3370,2902,2121,
    2635,3150,3920,3866,3690,3542,3430,3278];

AndersonDarling(BirthWeights, xlabel="Birth Weights");
# AD = 1.7168461277690312,  AD* = 1.748105851944578,  p value = 0.0001787724033071264
savefig("fig2.png")

R での結果

R"""
library(nortest)
options(digits=15)
a = ad.test($BirthWeights)
cat(a$statistic, a$p.value, "\n")
"""
# 1.71684612776903 0.000178772403307127 

検定のみを行うこともできる。

z = [76, 85, 93, 98, 100, 101, 111, 112, 119, 123];
AndersonDarling(z, plot=false);
# AD = 0.1671445157862479,  AD* = 0.18344110607540706,  p value = 0.9103125189292108

R"""
library(nortest)
options(digits=15)
a = ad.test($z)
cat(a$statistic, a$p.value, "\n")
"""
# 0.167144515786248 0.910312518929211