Julia に翻訳--010 対数軸でのプロット

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

対数軸でのプロット
http://aoki2.si.gunma-u.ac.jp/R/log-plot.html

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

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

takelog で "xy", "x", "y", "" を指定することにした。
デフォルトでは，普通の x-y グラフを描くことになる。

==========#

using Plots

function logplot(x, y; takelog = "", xlab = "x", ylab = "y", main = "", color = "gray")
    function logaxis(z)
        pos = zeros(100)
        minval = minimum(z)
        maxval = maximum(z)
        logmaxval = log10(maxval)
        logminval = log10(minval)
        width = logmaxval - logminval
        factor = logmaxval < 3.5 ? 15 : 8
        w = width / factor
        st = floor(logminval)
        delta = 10^st
        ic = 0
        while true
            for i = 1:9
                v = i * delta
                if minval - delta <= v <= maxval + delta
                    if (i == 1 || ic == 0) || (ic > 0 && log10(v) - pos[ic] > w)
                        if i == 1 && ic > 0 && log10(v) - pos[ic] <= w
                            ic -= 1
                        end
                        ic += 1
                        pos[ic] = log10(v)
                        prev = v
                    end
                end
            end
            delta *= 10
            if delta > maxval
                break
            end
        end
        return ic > 3 ? pos[2:ic-1] : [logminval, median([logminval, logmaxval]), logmaxval]
    end
    pyplot()
    logy = log10.(y)
    posy = logaxis(y)
    logx = log10.(x)
    posx = logaxis(x)
    if 'x' in takelog && 'y' in takelog
        plt = scatter(logx, logy, label = "")
        yticks!(posy, string.(round.(10 .^ posy, digits = 1)))
        xticks!(posx, string.(round.(10 .^ posx, digits = 1)))
    elseif 'y' in takelog
        plt = scatter(x, logy, label = "")
        yticks!(posy, string.(round.(10 .^ posy, digits = 1)))
    elseif 'x' in takelog
        plt = scatter(logx, y, label = "")
        xticks!(posx, string.(round.(10 .^ posx, digits = 1)))
    else
        plt = scatter(x, y, label = "")
    end
    plt = plot!(xlabel = xlab, ylabel = ylab, title = main)
    display(plt)
end

x = collect(1:0.5:9.5);
y = exp.(x);
logplot(x, y, main = "simple x : y graph")

logplot(x, y, takelog = "y", main = "semilog graph  x : log(y)")

logplot(y, x, takelog = "x", main = "semilog graph  log(x) : y")

x = 2 .^ collect(range(1, 5, length = 10));
y = 1.3 .^ collect(range(2, 10, length = 10));
logplot(x, y, takelog = "xy", main = "log-log graph  log(x) : log(y)")

Julia に翻訳--009 度数分布表を与えて，正規分布へのあてはめを行い，指定によっては図を描く

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

度数分布表を与えて，正規分布へのあてはめを行い，指定によっては図を描く
http://aoki2.si.gunma-u.ac.jp/R/fit_normal.html

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

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

図を描いた後に stdout に何らかの出力を追加すると，図が自動的には表示されない。
そのため，display() を億必要があるが，atom で実行していると図が 2 回出力されて煩わしい。

using で指定するパッケージに含まれる変数を，自分のプログラムでも使ってしまうと面倒なことになる。
==========#

using Distributions, Printf, Plots

function fitnormal(freq, lowerbound, width; accuracy = 0, method = "density",
                   xlab = "\$x\$", ylab = "\$f(x)\$", plotflag = true,
                   col = "cornsilk", col1 = "blueviolet", col2 = "chartreuse4",
                   verbose = true,)
    # method: "density" / "area"
    freq = vcat(0, freq, 0)
    m = length(freq)
    x = (lowerbound - accuracy / 2 - 3 * width / 2) .+ width .* collect(1:m)
    n = sum(freq)
    samplemean = sum(freq .* x) / n
    sd = sqrt(sum(freq .* x .^ 2) / n - samplemean^2)
    if method == "area"
        z = (x .+ (width / 2 - samplemean)) ./ sd
        F = cdf.(Normal(), z)
        F[end] = 1
        p = F .- vcat(0, F[1:end-1])
    else
        z = (x .- samplemean) ./ sd
        p = pdf.(Normal(), z) .* (width / sd)
    end
    expectation = n * p
    if verbose
        println("正規分布へのあてはめ (method = \"$method\")\n")
        @printf("推定された統計量　　平均値 = %.7g  標準偏差 = %.7g\n\n", samplemean, sd)
        println("   級中心    頻度      z 値       密度      期待値")
        for i = 1:size(x, 1)
            @printf("%8.4f %6i %9.3f %9.5f %10.2f\n",
                    x[i], freq[i], z[i], p[i], expectation[i])
        end
        @printf("     合計 %6i %19.5f %10.2f\n", sum(freq), sum(p), sum(expectation))
    end
    if plotflag
        pyplot()
        xl = (lowerbound - width / 2) .+ width .* collect(0:m)
        plt = bar(xl, vcat(freq, 0) ./ n, bar_width = width, fillcolor = col,
                  grid = false, xlabel = xlab, ylabel = ylab, label = "")
        scatter!(x, p, markercolor = col2, markerstrokecolor = col2, label = "")
        x2 = range(x[1], x[m], length = 100)
        plot!(x2, pdf(Normal(samplemean, sd), x2) * width, linecolor = col1, label = "")
        hline!([0], linecolor = "black", label = "")
        display(plt)
    end
    return Dict(:method => method, :mean => samplemean, :sd => sd, :x => x,
                :freq => freq, :z => z, :p => p, :expectation => expectation)
end

freq = [4, 19, 86, 177, 105, 33, 2]
a = fitnormal(freq, 40, 5)

正規分布へのあてはめ (method = "density")

推定された統計量　　平均値 = 57.98122  標準偏差 = 5.133511

   級中心    頻度      z 値       密度      期待値
 37.5000      0    -3.990   0.00014       0.06
 42.5000      4    -3.016   0.00412       1.75
 47.5000     19    -2.042   0.04833      20.59
 52.5000     86    -1.068   0.21974      93.61
 57.5000    177    -0.094   0.38686     164.80
 62.5000    105     0.880   0.26376     112.36
 67.5000     33     1.854   0.06964      29.67
 72.5000      2     2.828   0.00712       3.03
 77.5000      0     3.802   0.00028       0.12
     合計    426             0.99999     426.00



a = fitnormal(freq, 40, 5, method = "area")

正規分布へのあてはめ (method = "area")

推定された統計量　　平均値 = 57.98122  標準偏差 = 5.133511

   級中心    頻度      z 値       密度      期待値
 37.5000      0    -3.503   0.00023       0.10
 42.5000      4    -2.529   0.00549       2.34
 47.5000     19    -1.555   0.05428      23.12
 52.5000     86    -0.581   0.22070      94.02
 57.5000    177     0.393   0.37223     158.57
 62.5000    105     1.367   0.26129     111.31
 67.5000     33     2.341   0.07616      32.45
 72.5000      2     3.315   0.00915       3.90
 77.5000      0     4.289   0.00046       0.20
     合計    426             1.00000     426.00