Julia に翻訳--089 ポアソン分布への適合度の検定

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

ポアソン分布への適合度の検定
http://aoki2.si.gunma-u.ac.jp/R/poissondist.html

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

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

==========#

using Rmath, Printf, Plots

function poissondist(d, x)
    k = length(d)
    length(x) == k || error("度数ベクトル階級値ベクトルの長さは同じでなければなりません。")
    minx = minimum(x)
    maxx = maximum(x)
    o = zeros(maxx - minx + 1)
    o[x .- minx .+ 1] = d
    x = minx:maxx
    k = length(o)
    n = sum(o)
    lambda = sum(o .* x) / n
    p = dpois.(x, lambda)
    p[1] = ppois(minx, lambda)
    p[k] = ppois(maxx-1, lambda, false)
    e = n * p
    chisq, df, pvalue = calcpvalue(o, e)
    println("chisq. = $chisq,  df = $df,  p value = $pvalue")
    println("適合度")
    @printf("%10s  %4s  %10s  %10s\n", "階級", "度数", "確率", "期待値")
    for i = 1:k
        @printf("%10.6g  %4d  %10.6g  %10.6g\n", x[i], o[i], p[i], e[i])
    end
    pyplot()
    plt = bar(x, o, grid=false, tick_direction=:out, color=:blue, alpha=0.2,
        bar_width=1, xlabel="x", ylabel="Frequency", label="")
    scatter!(x, e, color=:blue, shape=:+, markersize=8, label="")
    display(plt)
    Dict(:chisq => chisq, :df => df, :pvalue =>pvalue)
end

function calcpvalue(o0, e0)
    o = copy(o0)
    e = copy(e0)
    while e[1] < 1
        o[2] += o[1]
        e[2] += e[1]
        popfirst!(o)
        popfirst!(e)
    end
    while e[end] < 1
        o[end-1] += o[end]
        e[end-1] += e[end]
        pop!(o)
        pop!(e)
    end
    chisq = sum((o .- e) .^ 2 ./ e)
    df = length(o) - 2
    p = pchisq(chisq, df, false)
    return chisq, df, p
end

d = [27, 61, 77, 71, 54, 35, 20, 11, 6, 2, 1];
x = 0:10;
poissondist(d, x)

diff(d)
diff(extrema(d))
f = [1, 2, 8, 8, 10, 15, 12, 15, 12, 6, 3, 1, 2, 0, 1]
x = 1:15
poissondist(f, x)

# chisq. = 14.143749368122279,  df = 8,  p value = 0.07809473292660538
# 適合度
#       階級  度数        確率      期待値
#          0    27    0.050198     18.3223
#          1    61    0.150181     54.8162
#          2    77    0.224655      81.999
#          3    71    0.224039     81.7743
#          4    54    0.167569     61.1627
#          5    35    0.100266     36.5971
#          6    20   0.0499957     18.2484
#          7    11    0.021368     7.79932
#          8     6  0.00799105     2.91673
#          9     2  0.00265639    0.969581
#         10     1  0.00108047    0.394373

Julia に翻訳--088 正規分布への適合度の検定

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

正規分布への適合度の検定
http://aoki2.si.gunma-u.ac.jp/R/normaldist.html

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

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

関数の引数として受け取った配列（ベクトル）を変更すると，その結果が呼出元へ反映される。
それって，R とは違い，便利なときもあるし不便な（困る）ときもある。

==========#

using Rmath, Printf, Plots

function normaldist(x, b, w, a)
    n = sum(x)
    x = vcat(0, x, 0)
    k = length(x)
    mid = (b-w/2:w:b+k*w-w) .- a ./ 2
    xbar = sum(mid .* x) / n
    variance = sum(x .* (mid .- xbar) .^ 2) / n
    SD = sqrt(variance)
    println("平均値 = $xbar,  不偏分散 = $variance,  標準偏差 = $SD")
    z = ((mid .+ w/2) .- xbar) / SD
    p = pnorm.(z)
    p[k] = 1
    p = p .- vcat(0, p[1:end-1])
    expectation = n * p
    chisq, df, pvalue = calcpvalue(x, expectation)
    println("chisq. = $chisq,  df = $df,  p value = $pvalue")
    println("適合度")
    @printf("%10s  %4s  %10s  %10s  %10s\n", "級中心", "度数", "標準化得点", "確率", "期待値")
    for i = 1:k
        @printf("%10.6g  %4d  %10.6g  %10.6f  %10.6g\n",
                mid[i], x[i], z[i], p[i], expectation[i])
    end
    plt = bar(mid, x, grid=false, tick_direction=:out, color=:blue, alpha=0.2,
        bar_width=w, xlabel="x", ylabel="Frequency", label="")
    x3 = range(mid[1]-w, mid[k]+w, length=1000)
    y = dnorm.(x3, xbar, SD) .* (w * n)
    plot!(x3, y, color=:black, label="")
    y = dnorm.(mid, xbar, SD) .* (w * n)
    scatter!(mid, y, color=:blue, shape=:+, markersize=8, label="")
    display(plt)
    savefig(plt, "fig1.png")
    Dict(:chisq => chisq, :df => df, :pvalue =>p,
        :Mean => xbar, :Variance => variance, :SD => SD)
end

function calcpvalue(x, expectation)
    x2 = copy(x) # 絶対に必要
    e2 = copy(expectation) # 絶対に必要
    while e2[1] < 1
        x2[2] += x2[1]
        e2[2] += e2[1]
        popfirst!(x2)
        popfirst!(e2)
    end
    while e2[end] < 1
        x2[end-1] += x2[end]
        e2[end-1] += e2[end]
        pop!(x2)
        pop!(e2)
    end
    chisq = sum((x2 .- e2) .^ 2 ./ e2)
    df = length(x2) - 3
    pvalue = pchisq(chisq, df, false)
    return chisq, df, pvalue
end

x = [4, 19, 86, 177, 105, 33, 2]
res = normaldist(x, 40, 5, 0.1)

# 平均値 = 57.931220657276995,  不偏分散 = 26.352933721263412,  標準偏差 = 5.133510857226603
# chisq. = 6.000195029641498,  df = 4,  p value = 0.19913370901373778
# 適合度
#     級中心  度数  標準化得点        確率      期待値
# k = 9
#      37.45     0    -3.50271    0.000230   0.0980958
#      42.45     4    -2.52872    0.005494      2.3403
#      47.45    19    -1.55473    0.054281     23.1238
#      52.45    86   -0.580737    0.220704     94.0197
#      57.45   177    0.393255    0.372226     158.568
#      62.45   105     1.36725    0.261292      111.31
#      67.45    33     2.34124    0.076164     32.4459
#      72.45     2     3.31523    0.009152     3.89878
#      77.45     0     4.28922    0.000458    0.195038

Julia に翻訳--087 二項分布への適合度の検定

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

二項分布への適合度の検定
http://aoki2.si.gunma-u.ac.jp/R/binomdist.html

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

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

d = d[-1] なんてのを見落としそうになった。

==========#

using Rmath, Printf, Plots

function binomdist(d, x, size)
    k = length(d)
    k == length(x) || error("度数ベクトル d と階級値ベクトル x の長さが違います")
    all(floor.(Int, d) .== d)|| error("度数ベクトル中に小数値があります")
    all(d .>= 0) || error("度数ベクトル中に負の値があります")
    all(x .<= size) ||error("階級値ベクトル中に試行数より大きい数値があります")
    all(x .>= 0) || error("階級値ベクトル中に負の数値があります")
    all(floor.(Int, x) == x) || error("階級値ベクトル中に小数値があります")
    n = sum(d)
    prob = sum(d .* x) / n / size
    p = dbinom.(x, size, prob)
    e = n * p
    chisq, df, pvalue = calcpvalue(d, e)
    println("chisq. = $chisq,  df = $df,  p value = $pvalue")
    println("適合度")
    @printf("%4s  %4s  %10s  %10s\n", "階級", "度数", "確率", "期待値")
    for i = 1:k
        @printf("%4d  %4d  %10.6f  %10.4f\n", x[i], d[i], p[i], e[i])
    end
    pyplot()
    plt = bar(x, d, grid=false, tick_direction=:out, color=:blue, alpha=0.2,
         xlabel="x", ylabel="Frequency", label="")
    scatter!(x, e, color=:blue, shape=:+, markersize=8, label="")
    display(plt)
    Dict(:chisq => chisq, :df => df, :pvalue => pvalue, :probability => prob)
end

function calcpvalue(d0, e0)
    d = copy(d0)
    e = copy(e0)
    while e[1] < 1
        d[2] += d[1]
        e[2] += e[1]
        popfirst!(d)
        popfirst!(e)
    end
    while e[end] < 1
        d[end-1] += d[end]
        e[end-1] += e[end]
        pop!(d)
        pop!(e)
    end
    chisq = sum((d .- e) .^ 2 ./ e)
    df = length(d) - 2
    pvalue = pchisq(chisq, df, false)
    return chisq, df, pvalue
end

d = [2, 14, 20, 34, 22, 8];
x = 0:5
size = 5
binomdist(d, x, 5)
# chisq. = 4.007575353728885,  df = 4,  p value = 0.4049816078597148
# 適合度
# 階級  度数        確率      期待値
#    0     2    0.015046      1.5046
#    1    14    0.098913      9.8913
#    2    20    0.260105     26.0105
#    3    34    0.341989     34.1989
#    4    22    0.224826     22.4826
#    5     8    0.059121      5.9121