#==========
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)
