算額あれこれ

算額問題をコンピュータで解きます

Julia に翻訳--006 正規確率紙に累積相対度数をプロットする

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

正規確率紙に累積相対度数をプロットする
http://aoki2.si.gunma-u.ac.jp/R/npp2.html

npp2 <- function(x) # データベクトル
{
        x <- x[!is.na(x)] # 欠損値を持つケースを除く
        n <- length(x) # データの個数
        x <- sort(x) # 昇順に並べ替える
        y <- (1:n-0.5)/n # 累積相対度数
        probs <- c(0.01, 0.1, 1, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 99, 99.9, 99.99)/100
        plot(c(x[1], x[n]), qnorm(c(probs[1], probs[17])), type="n", yaxt="n",
                xlab="Observed Value", ylab="Cumulative Percent",
                main="Normal Probability Paper")
        points(x, qnorm(y))
        axis(2, qnorm(probs), probs*100)
}

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

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

npp() よりは簡単?
指数乱数は randexp() だった。R はその他の分布型の乱数も用意されているのになあ。

==========#

using Distributions, Plots, Plots.PlotMeasures

function npp2(x; leftmargin=2mm, xlab = "観察値", ylab = "累積確率",
        main = "正規確率紙")
    pyplot()
    x = collect(skipmissing(x))
    n = length(x)
    sort!(x)
    y = (collect(1:n) .- 0.5) ./ n;
    qny = quantile(Normal(), y);
    probability = [0.01, 0.1, 1, 5, 10, 20, 30, 40, 50,
        60, 70, 80, 90, 95, 99, 99.9, 99.99]/100;
    qnp = quantile(Normal(), probability)
    scatter(x, qny,
        tickdirection=:out,
        grid=false,
        label="",
        xlabel=xlab,
        ylabel=ylab,
        title=main,
        leftmargin=leftmargin,
        ylims=(minimum(qny), maximum(qny)),
        size=(600, 400))
    yticks!(qnp, string.(round.(probability, digits=4)))
    hline!(qnp, linecolor = :gray, linestyle=:dot, label="")
end

using Random

Random.seed!(123)
x = randn(100) # 正規乱数
npp2(x)

Random.seed!(123)
npp2(rand(100)) # 一様乱数

Random.seed!(123)
npp2(randexp(100)) # 指数乱数