Julia に翻訳--008 測定値をワイブル確率紙にプロットする

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

測定値をワイブル確率紙にプロットする
http://aoki2.si.gunma-u.ac.jp/R/weibull.html

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

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

対数正規確率紙を書いた後だから，簡単だった。

==========#

using Distributions, Plots, Plots.PlotMeasures

function weibull(x; leftmargin = 2mm, xlab = "観察値", ylab = "累積確率", main = "ワイブル確率紙")
    weib(p) = log10(log10(1 / (1 - p)))
    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()
    n = length(x)
    logx = log10.(sort!(x))
    y = weib.((collect(1:n) .- 0.5) ./ n)
    y0 = vcat(
        10.0 .^ (-10:0),
        [5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 99, 99.9, 99.99, 99.999])
    y0 = y0[y0.>10/n]
    probability = weib.(y0 / 100)
    scatter(logx, y, label = "", xlabel = xlab, ylabel = ylab, title = main)
    yticks!(probability, string.(round.(y0, digits = 4)))
    pos = logaxis(x)
    xticks!(pos, string.(round.(10 .^ pos, digits = 1)))
end

x = [
    1.80531046285254,3.19311612835466,3.28834436413816,6.09849913246918,1.38753827431415,
    4.25956917997002,2.19512553656654,0.512283970489512,4.34705764607181,1.11953153395821,
    7.68555210540847,1.7643219025482,0.986773337068037,0.682127148920385,2.6323987051625,
    3.08191138366997,0.575444181661929,7.04322086787495,0.867594222227101,1.24874964755851,
    0.618424866491301,3.98775421377481,0.601446810729256,3.45839946486477,1.52155116044975,
    6.673646245403,1.1408427257711,0.231056233767167,4.62036593935406,0.568213427986162,
    2.49540860834656,1.18903624694958,0.886439743171324,2.1480855169868,4.52228064341519,
    1.15375714297612,5.86131383827675,0.215763115606156,3.32187592926364,0.542381907439051,
    2.47025306148576,2.40605148698434,2.84101535087403,1.52931862347578,0.669626351294388,
    3.40654318414194,1.77750271751872,5.56339366874868,0.626310926516198,6.24842979611537,
    3.15329814448187,2.97783388594968,4.72055464574222,4.97361132035854,2.71284323516702,
    3.84714124059119,0.411863928700296,1.25628023322285,1.43876907120354,2.90263414177699,
    1.52999869396695,0.933096212140327,2.62990034345694,0.765430620280283,3.39133399716817,
    1.80780059565393,1.37695181832508,2.64438749146156,1.21643677550211,1.72673120016061,
    3.31002652328123,1.00147436267977,3.12131155688389,3.17198635435594,2.21115310778158,
    1.15469936237419,2.50724717782027,4.06120539534992,5.13987643842686,4.78868160713468,
    0.837216840665425,0.755663424463497,4.18810116764077,6.56524669360705,6.42437126634158,
    6.28044699738115,1.08706640499649,1.84654965091899,2.3296569581759,2.90964806098607,];
weibull(x)

using Random
Random.seed!(123);
x = rand(20) .* 100 .+ 2.1;
weibull(x)

Julia の小ネタ--005 : は ^ より演算順位が低い

10^-10, 10^-9, ..., 1 の数列を得ようとして

10 .^ -10:0 としても，1.0000000000000006e-10:1.0:-0.9999999999 になるだけ。これは何も含まない。
float.(1.0000000000000006e-10:1.0:-0.9999999999) は Float64[] で「空」である。: は ^ より演算順位が低いからである（R は逆）

それならば，これはどうだ。

10 .^ (-10:0) はエラーになる。

DomainError with -10:
Cannot raise an integer x to a negative power -10.
Make x or -10 a float by adding a zero decimal (e.g., 2.0^-10 or 2^-10.0 instead of 2^-10), or write 1/x^10, float(x)^-10, x^float(-10) or (x//1)^-10

解決策のサジェスチョンどおり，以下のようにすると

10.0 .^ (-10:0) または 10 .^ (-10.0:0) でやっと望みの結果が得られた。やれやれ。

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

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

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

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

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

軸を対数目盛りにする :log はあるが，表示が 10^n などとなるので，指数表示ではなく元の数値の表示に従うように新たに軸メモリを計算する関数を作ったが，これが手こずった。

==========#

using Distributions, Plots, Plots.PlotMeasures

function lnpp(x; leftmargin = 2mm, xlab = "観察値", ylab = "累積確率",
              main = "対数正規確率紙")
    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()
    n = length(x)
    logx = log10.(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)
    margin = (log10(maximum(x)) - log10(minimum(x))) / 600 * 30
    scatter(
        logx,
        qny,
        label = "",
        xlims = (log10(minimum(x)) - margin, log10(maximum(x)) + margin),
        ylims = (minimum(qny), maximum(qny)),
        xlabel = xlab,
        ylabel = ylab,
        title = main,
    )
    yticks!(qnp, string.(round.(probability, digits = 4)))
    hline!(qnp, linecolor = :gray, linestyle = :dot, label = "")
    pos = logaxis(x)
    xticks!(pos, string.(round.(10 .^ pos, digits = 1)))
end

using Random
Random.seed!(123)
x = randn(100) .* 10 .+ 80
lnpp(x)

x = [
    32.4978550221931,3.11052584094339,9.55400310478408,30.3645286610177,
    21.1624392829688,7.37476464691322,28.9129834585325,1.8042322428236,
    8.99964598570283,3.55895261576496,9.62585749164628,6.35726278028454,
    1.95005492569664,13.4876154673338,8.24295995200768,6.48532599314323,
    4.35588525758088,19.7888507959163,2.54536615538532,22.1792645172972,
    22.6290782281654,11.394279463117,1.62423013916794,1.61310795035265,
    32.1156725493402,24.4324047816443,4.62558675355759,5.9662576846906,
    1.67772428988997,3.96909669181296,12.3842951428101,4.75229826175545,
    3.14506115407662,10.8612899803449,3.24849536568508,11.4524807755855,
    8.44567923387873,7.53737447514118,16.0361633159583,30.5856770679345,
    34.2904590210926,37.4844432670513,25.5460725535086,5.04602127463469,
    3.10504378339411,28.8353322772864,7.37368663949144,42.5958262784602,
    0.946681530550752,9.42924224403108,]
lnpp(x)