Julia に翻訳--034 AIC による，ヒストグラム（度数分布表）の最適階級分割の探索

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

AIC による，ヒストグラム（度数分布表）の最適階級分割の探索
http://aoki2.si.gunma-u.ac.jp/R/AIC-Histogram.html

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

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

ほぼそのまま書き写すことができた。

==========#

using Statistics, DataFrames, Plots

function aichistogram(x; d = 0, c = floor(Int, 2*sqrt(length(x))-1))
    function rect(x1, y1, x2, y2; col=:black, lty=:solid, lwd=1, fcol="", alpha=0.2)
        if fcol == ""
            plot!([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], linecolor=col,
                  linestyle=lty, linewidth=lwd, label="")
        else
            plot!([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], linecolor=col,
                  linestyle=lty, seriestype=:shape, linewidth=lwd, fillcolor=fcol,
                  label="", alpha=alpha)
        end
    end

    function model(c1, r, c2, n, c)
        logNZ(x, y) = x > 0 ? x*log(x/y) : 0
        N = sum(n)
        sum1 = sum(n[1:c1])
        sum2 = sum(n[(c-c2+1):c])
        cc1c2r1 = (c-c1-c2) ÷ r + 1
        temp = 0
        for j in 2:cc1c2r1
            sum3 = sum(n[(c1+(j-2)*r+1):(c1+(j-1)*r)])
            temp = temp+logNZ(sum3, r*N)
        end
        AIC = -2*(logNZ(sum1, c1*N)+temp+logNZ(sum2, c2*N))+2*cc1c2r1
        return(AIC)
    end

    N = length(x)
    minx, maxx = extrema(x)
    brks = range(minx-d/2, maxx+d/2, length=c+1)
    lo = brks[1]
    hi = brks[end]
    width = brks[2] - brks[1]
    x2 = floor.(Int, (x .- minx) ./ width) .+ 1
    maxx = maximum(x2)
    n = zeros(Int, maxx)
    for i = 1:N
        n[x2[i]] += 1
    end
    df = DataFrame(c1=0 , r=0, c2=0, AIC=Inf)
    for c1 = 1:c-2
        for c2 = 1:c-2
            if c <= c1+c2
                continue
            end
            for r in 1:c-c1-c2
                if (c-c1-c2) % r == 0
                    append!(df,
                      DataFrame(c1=c1, r=r, c2=c2, AIC=model(c1, r, c2, n, c)))
                end
            end
        end
    end
    p = n / N * 100
    sort!(df, 4)
    pyplot()
    c1  = df[1, 1]
    r   = df[1, 2]
    c2  = df[1, 3]
    AIC = df[1,4]
    plt = bar(brks, n, grid=false, bar_width=width, color=:blue, alpha=0.1,
                label="")
    m = floor(Int, (c-c1-c2)/r)
    cmg = cumsum(vcat(0, c1, repeat([r], m), c2))
    for i = 1:m+2
        rect(lo+cmg[i]*width, 0,
             lo+cmg[i+1]*width, mean(p[(cmg[i]+1):(cmg[i+1])]),
             fcol=:red, alpha=0.5)
    end
    annotate!(hi, maximum(n),
              text("AIC =" * string(round(AIC, digits=2)), :top, :right, 10))
    display(plt)
    Dict(:df => df, :n => n, :breaks => brks)
end

x = [28.67, 40.29, 10.61, 33.85, 36.19, 20.63, 9.64, 15.26,
  15.53, 73.62, 63.29, 32.77, 32.28, 11.90, 54.16, 4.73, 24.67,
  17.66, 25.84, 22.89, 15.68, 5.48, 36.41, 20.33, 44.58, 57.23,
  65.89, 57.91, 2.39, 9.15, 10.27, 3.04, 12.35, 32.78, 44.23,
  31.14, 6.03, 27.90, 28.73, 42.09, 3.99, 9.74, 6.85, 0.16, 9.26,
  7.72, 34.42, 32.77, 6.80, 10.45, 29.80, 5.89, 13.56, 50.55, 0.51,
  0.19, 7.19, 5.94, 11.24, 32.32, 15.27, 29.64, 10.03, 2.01, 13.89,
  20.83, 27.49, 14.46, 8.22, 27.81, 33.65, 38.57, 8.66, 1.40,
  23.97, 15.11, 63.32, 7.76, 1.58, 48.66, 44.46, 0.02, 38.12,
  18.51, 101.75, 34.16, 27.99, 5.22, 1.82, 8.22, 4.89, 97.50, 2.10,
  26.19, 10.11, 8.39, 25.83, 1.05, 25.63, 18.35]

  a = aichistogram(x, d = 0.01)

Julia に翻訳--033 度数分布表

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

度数分布表
http://aoki2.si.gunma-u.ac.jp/R/dosuu-bunpu.html

ファイル名: dosuubunpu.jl　　関数名: dosuubunpu 同じ名前の 3 関数

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

度数分布表作成は FreqTables の freqtable を使うのが便利だが，一応書いてみた。

==========#

using Printf, Plots

function printtable(category, freq, pcnt, cumpcnt)
    @printf("%10s %7s %11s  %11s\n", "階級", "度数", "パーセント", "累積パーセント")
    for i = 1:length(category)
        @printf("%10s %7d %11.1f  %11.1f\n", category[i], freq[i], pcnt[i], cumpcnt[i])
    end
end

function db(category, freq, percent, color, alpha, title, ylabel, xlabel)
    pcnt = freq ./ (100.0 / sum(freq))
    cumpcnt = cumsum(pcnt)
    pyplot()
    if percent
        yvalue = copy(pcnt)
        ylabel == "Frequency" && (ylabel = "Percent")
    else
        yvalue = copy(freq)
    end
    plt = bar(category, yvalue, tick_direction=:out, grid=false, bar_width=1,
                 linewidth=0.5, color=color, alpha=alpha, title=title,
                ylabel=ylabel, xlabel=xlabel, label="")
    display(plt)
    printtable(string.(category), freq, pcnt, cumpcnt)
    Dict(:category => category, :freq => freq, :pcnt => pcnt, :cumpcnt => cumpcnt)
end

# 離散分布データ（整数型）の場合（度数 0 も含む度数分布表を作成し，棒グラフを描く）
function dosuubunpu(x::Vector{Int64}; percent::Bool = false,
                    color=:azure2, alpha=0.5,
                    xlabel="", ylabel="Frequency", title="")
    n = length(x)
    mn, mx = extrema(x)
    len = mx - mn + 1
    freq = zeros(Int, len)
    for i = 1:n
        freq[x[i] - mn + 1] += 1
    end
    category = collect(mn:mx)
    db(category, freq, percent, color, alpha, title, ylabel, xlabel)
end

# 離散分布データ（文字列型）の場合（度数分布表を作成し，棒グラフを描く）
function dosuubunpu(x::Vector{String}; percent::Bool = false, color=:honeydew,
                    alpha=0.5, xlabel="", ylabel="Frequency", title="")
    category = [s for s in Set(x)]
    sort!(category)
    freq = zeros(Int, length(category))
    for (i, cat) in enumerate(category)
        freq[i] = count(isequal(cat), x)
    end
    db(category, freq, percent, color, alpha, title, ylabel, xlabel)
end

# 連続分布データの場合（ヒストグラムを描く）
function dosuubunpu(x::Vector{Float64}, min::Float64, w::Float64;
                    color=:floralwhite, ylabel="", density = false, title="")
    n = length(x)
    mn, mx = extrema(x)
    while mn < min
        min -= w
    end
    len = floor(Int, (mx - min) / w)
    freq = zeros(len)
    for i =1:n
        freq[floor(Int, (x[i] - min)/w)] += 1
    end
    category = collect(1:len) .* w .+ min
    pcnt = freq ./ n .* 100.0
    cumpcnt = cumsum(pcnt)
    if ylabel == "" && density
        normalize = :pdf; ylabel = "Density"
    elseif ylabel == ""
        normalize = false; ylabel = "Frequency"
    end
    pyplot()
    plt = histogram(x, bins=vcat(category, category[end]+w),
        fillcolor = color, linewidth = 0.5, tick_direction = :out,
        normalize = normalize, ylabel = ylabel, title=title, label="")
    display(plt)
    printtable(string.(category), freq, pcnt, cumpcnt)
    Dict(:category => category, :freq => freq, :pcnt => pcnt, :cumpcnt => cumpcnt)
end

x = [
    162, 159, 163, 157, 152, 168, 153, 156, 167, 161,
    154, 162, 160, 157, 169, 160, 162, 158, 161, 160,
    163, 160, 163, 153, 164, 163, 163, 153, 155, 155,
    162, 163, 168, 160, 158, 168, 163, 163, 158, 153,
    161, 153, 168, 156, 155, 159, 158, 161, 157, 155,
    161, 156, 167, 156, 158, 152, 160, 160, 155, 157,
    158, 160, 157, 156, 164, 157, 161, 158, 161, 153,
    163, 161, 160, 162, 159, 162, 161, 158, 160, 173
    ];

データが実数型の場合はヒストグラム
a = dosuubunpu(float.(x), 150.0, 2.0);
b = dosuubunpu(float.(x), 150.0, 2.0, density=true);

整数型の場合には棒グラフ
c = dosuubunpu(x)
c = dosuubunpu(x, percent=true)

文字列型の場合も棒グラフ
y = [
    "k", "h", "l", "f", "a", "o", "b", "e", "n", "j",
    "c", "k", "i", "f", "p", "i", "k", "g", "j", "i",
    "l", "i", "l", "b", "m", "l", "l", "b", "d", "d",
    "k", "l", "o", "i", "g", "o", "l", "l", "g", "b",
    "j", "b", "o", "e", "d", "h", "g", "j", "f", "d",
    "j", "e", "n", "e", "g", "a", "i", "i", "d", "f",
    "g", "i", "f", "e", "m", "f", "j", "g", "j", "b",
    "l", "j", "i", "k", "h", "k", "j", "g", "i", "q"
    ];

d = dosuubunpu(y)
d = dosuubunpu(y, percent=true)