Julia に翻訳--029 ワイブル分布のパラメータ（最尤推定）

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

ワイブル分布のパラメータ（最尤推定）
http://aoki2.si.gunma-u.ac.jp/R/weibull-par.html

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

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

特に問題なし。

==========#

function weibullpar(x; epsilon = 1e-7, maxloop = 500)
    n = length(x)
    m = a = 1
    for i = 1:maxloop
        ao = a
        mo = m
        temp1 = x .^ mo
        temp2 = log.(x)
        a = n / sum(temp1)
        m = n / (a * sum(temp1 .* temp2) - sum(temp2))
        if abs(a - ao) < epsilon && abs(m - mo) < epsilon
            scale = (1 / a) ^ (1 / m)
            return Dict(:m => m, :scale => scale)
        end
    end
    error("not comberged")
end

x = [
    2.9614308, 2.7978909, 1.2497852, 3.2014155, 2.1161621, 0.4623305, 
    2.9601025, 1.40131, 7.5831998, 3.7555033, 5.0130449, 1.1536598, 
    2.7115764, 1.6891542, 4.5575992, 1.062815, 2.5742121, 3.6710305, 
    1.7050176, 1.9851916
    ];
weibullpar(x) # (1.8033671659326795, 3.083428778727086)
# Dict{Symbol,Float64} with 2 entries:
#   :m     => 1.80337
#   :scale => 3.08343

 

Julia に翻訳--028 ポリヤ・エッゲンベルガー分布

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

ポリヤ・エッゲンベルガー分布
http://aoki2.si.gunma-u.ac.jp/R/PolyaEggenberger.html

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

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

引数の個数が違うので，同じ名前で 2 種類の関数定義ができる。

==========#

using SpecialFunctions # logfactorial

function polyaeggenberger(x::Int64, n::Int64, p::Float64, δ::Float64)
    logcomb(n, i) = logfactorial(n) - logfactorial(i) - logfactorial(n - i)
    exp(
        logcomb(n, x) +
        sum([i == 0 ? 0 : log(p + (i - 1) * δ) - log(1 + (i - 1) * δ) for i = 0:x]) +
        sum([log(1 - p + (i - x - 1) * δ) - log(1 + (i - 1) * δ) for i = x+1:n])
    )
end

# λ, r を与える場合
function polyaeggenberger(x::Int64, λ::Float64, r::Float64)
    exp(                                            # 対数で計算して最後に逆対数をとる
        sum([i < 0 ? 0 : log(1 + i * r / λ) for i = 0:x-1]) +
        x * log(λ) - logfactorial(x) +
        (-λ / r - x) * log(1 + r)
    )
end

n = 800 # n, p, δ を与える場合
p = 0.00373625
δ = 0.000322
polyaeggenberger(0, n, p, δ) # 0.06964383745279974
polyaeggenberger(1, n, p, δ) # 0.16606182454351842
polyaeggenberger(2, n, p, δ) # 0.21483159643906802
polyaeggenberger(3, n, p, δ) # 0.19978508436177314
sum(polyaeggenberger(x, n, p, δ) for x in 0:15) # 0.9999914961492099

λ = 2.989 # λ, r を与える場合
r = 0.2576
polyaeggenberger(0, λ, r) # 0.06998131064894278
polyaeggenberger(1, λ, r) # 0.1663280355675015
polyaeggenberger(2, λ, r) # 0.21469489514688703
polyaeggenberger(3, λ, r) # 0.1994099479362071
sum(polyaeggenberger(x, λ, r) for x = 0:15) # 0.9999907914029964

Julia に翻訳--027 負の超幾何分布

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

負の超幾何分布
http://aoki2.si.gunma-u.ac.jp/R/negative-geometric.html

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

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

ガンマ関数の自然対数を使って階乗を計算
ガンマ関数の代わりに logfactorial を使おう

==========#

using SpecialFunctions # logfactorial

function negativegeometric(x, N, n, r)
    logcomb(n, i) = logfactorial(n) - logfactorial(i) - logfactorial(n - i)
    exp(logcomb(x - 1, r - 1) + logcomb(N - x, n - r) - logcomb(N, n))
end

negativegeometric(4, 5000, 400, 3) # 0.0014042237146677339

Julia に翻訳--026 多項分布

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

多項分布
http://aoki2.si.gunma-u.ac.jp/R/multinomial.html

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

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

logfactorial(n) ≡ loggamma(n + 1)

==========#

using SpecialFunctions # logfactorial

function multinomial(x::Array{Int64,1}, p::Array{Float64,1})
    length(x) == length(p) || error("length(x) != length(p)")
    sum(p) == 1            || error("sum(p) != 1")
    all(x .>= 0)           || error("any(x) < 0")
    exp(logfactorial(sum(x)) + sum(x .* log.(p)) - sum(logfactorial.(x)))
end

x = [4, 2, 3, 0];
p = [0.5, 0.15, 0.3, 0.05];
multinomial(x, p) # 0.047840625000000046

Julia にもある。使い方が他の統計関数と同じく Julia 独特？

using Distributions
n = sum(x)
pdf(Multinomial(n, p), [4, 2, 3, 0]) # 0.047840625000000025

 

Julia に翻訳--025 プロット関数群

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

プロット関数群
http://aoki2.si.gunma-u.ac.jp/R/plot.html

ファイル名: plotter.jl　　関数名: plot*

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

なかなか楽しかった。色々拡張した。
前に NHK でやっていた，家紋を描くシリーズを見ていたので，ちょっとやってみた。

==========#

using Plots

# 開始
function plotbegin(; w = 600, h = 600)
    pyplot()
    plt = plot(xlabel="", ylabel="", grid=false, showaxis=false,
               aspect_ratio=1, size=(w, h), label="")
end

# 終了
function plotend()
    plot!()
end

# （x1, y1) - (x2, y2) の線分
function plotline(x1, y1, x2, y2; col=:black, lty=:solid, lwd=1)
    plot!([x1, x2], [y1, y2], linecolor=col, linestyle=lty, linewidth=lwd, label="")
end

# （x1, y) - (x2, y) の水平線分
function plothline(x1, x2, y; col=:black, lty=:solid, lwd=1)
    plot!([x1, x2], [y, y], linecolor=col, linestyle=lty, linewidth=lwd, label="")
end

# (x, y1) - (x, y2) の垂直線分
function plotvline(x, y1, y2; col=:black, lty=:solid, lwd=1)
    plot!([x, x], [y1, y2], linecolor=col, linestyle=lty, linewidth=lwd, label="")
end

# (x1, y1), (x2, y2) を対角座標とする矩形
function plotbox(x1, y1, x2, y2; col=:black, lty=:solid, lwd=1, fcol="")
    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="")
    end
end

# (ox, oy) を原点とする半径 r の円（円弧）。fcol を指定すると塗りつぶし。
function plotcircle(ox, oy, r; startangle=0, endangle=360,
                    col=:black, lty=:solid, lwd=1, fcol="")
    plotellipse(ox, oy, r, r, φ=0, startangle=startangle,
                endangle=endangle, col=col, lty=lty, lwd=lwd, fcol=fcol)
end

# 度をラジアンに変換する
radian(degree) = degree / 180 * π

# (ox, oy) を中心，ra, rb を半径とする楕円（楕円弧）。fcol を指定すると塗りつぶし。
function plotellipse(ox, oy, ra, rb; φ=0, startangle=0, endangle=360,
                     length=100, col=:black, lty=:solid, lwd=1, fcol="")
    θ = vcat(range(radian(startangle), radian(endangle),
              length=length), radian(endangle))
    if φ == 0
        if fcol == ""
            plot!(ra .* cos.(θ) .+ ox, rb .* sin.(θ) .+ oy,
                  linecolor=col, linestyle=lty, linewidth=lwd, label="")
        else
            plot!(ra .* cos.(θ) .+ ox, rb .* sin.(θ) .+ oy,
                  linecolor=col, linestyle=lty, linewidth=lwd,
                  seriestype=:shape, fillcolor=fcol, label="")
        end
    else
        x = ra .* cos.(θ)
        y = rb .* sin.(θ)
        φ = radian(φ)
        cosine = cos(φ)
        sine = sin(φ)
        if fcol == ""
            plot!(cosine .* x .- sine .* y .+ ox,
                  sine .* x .+ cosine .* y .+ oy,
                  linecolor=col, linestyle=lty, linewidth=lwd, label="")
        else
            plot!(cosine .* x .- sine .* y .+ ox,
                  sine .* x .+ cosine .* y .+ oy,
                  linecolor=col, linestyle=lty, linewidth=lwd,
                  seriestype=:shape, fillcolor=fcol, label="")
        end
    end
end

# 任意の多角形。fcol を指定すると塗りつぶし。
function plotpolygon(x, y; col=:black, lty=:solid, lwd=1, fcol="")
    x = vcat(x, x[1])
    y = vcat(y, y[1])
    if fcol == ""
        plot!(x, y, linecolor=col, linestyle=lty, linewidth=lwd, label="")
    else
        plot!(x, y, linecolor=col, linestyle=lty, linewidth=lwd,
              seriestype=:shape, fillcolor=fcol, label="")
    end
end

# (x1, y1) を描き始め，一辺の長さ l の正 n 多角形。fcol を指定すると塗りつぶし。
function plotpolygon1(x1, y1, l, n; col=:black, lty=:solid, lwd=1, fcol="")
    θ = range(0, 2π, length=n+1)
    x = fill(float(x1), n+1)
    y = fill(float(y1), n+1)
    for i = 2:n
        x[i] = x[i-1]+l*cos(θ[i])
        y[i] = y[i-1]+l*sin(θ[i])
    end
    if fcol == ""
        plot!(x, y, linecolor=col, linestyle=lty, linewidth=lwd, label="")
    else
        plot!(x, y, linecolor=col, linestyle=lty,
              seriestype=:shape, fillcolor=fcol, linewidth=lwd, label="")
    end
end

# (ox, oy) を中心，r を半径とする円に内接する正 n 多角形。fcol を指定すると塗りつぶし。
function plotpolygon2(ox, oy, r, n; φ=90, col=:black, lty=:solid, lwd=1, fcol="")
    θ = range(0, 2π, length=n+1) .+ radian(φ)
    if fcol == ""
        plot!(r .* cos.(θ) .+ ox, r .* sin.(θ) .+ oy,
              linecolor=col, linestyle=lty, linewidth=lwd, label="")
    else
        plot!(r .* cos.(θ) .+ ox, r .* sin.(θ) .+ oy,
              linecolor=col, linestyle=lty, linewidth=lwd,
              seriestype=:shape, fillcolor=fcol, label="")
    end
end

# (x1,y1) - (x2, y2) の矩形範囲に格子を描く。間隔は wx, wy。
function plotgrid(x1, y1, x2, y2, wx; wy=wx, col=:gray, lty=:dot, lwd=1)
    X1 = min(x1, x2)
    X2 = max(x1, x2)
    Y1 = min(y1, y2)
    Y2 = max(y1, y2)
    for i = 0:fld(abs(X2-X1), wx)
        plot!([X1+wx*i, X1+wx*i], [Y1, Y2], linecolor=col,
              linestyle=lty, linewidth=lwd, label="")
    end
    for i = 0:fld(abs(y2-y1), wy)
        plot!([X1, X2], [Y1+wy*i, Y1+wy*i], linecolor=col,
              linestyle=lty, linewidth=lwd, label="")
    end
end

# ネフロイド
plotbegin()
t = 0:3:360
l = 130
ox = l .* cos.(radian.(t))
oy = l .* sin.(radian.(t))
for i = 1:length(t)
    plotcircle(ox[i]+250, oy[i]+250, abs(ox[i]),
               startangle=0, endangle=360, col="blue")
end
plotend()

# カージオイド
plotbegin()
t = 0:5:360
l = 100
ox = l .* cos.(radian.(t))
oy = l .* sin.(radian.(t))
for i = 1:length(t)
    plotcircle(ox[i]+320, oy[i]+250,
               sqrt((ox[i]-ox[1])^2+(oy[i]-oy[1])^2),
               startangle=0, endangle=360, col="blue")
end
plotend()

# 多角形1
using FixedPointNumbers, ColorTypes # for RGB
plotbegin();
for i = 30:-1:3
    plotpolygon1(0, 0, 70, i, fcol=RGB((70-i)/70, 0, 0))
end
plotend()

# 多角形2
plotbegin();
for i = 30:-1:3
    plotpolygon2(0, 0, 70, i, fcol=RGB(0, (30-i)/30, 0))
end
plotend()

plotbegin(h=300)
plotpolygon2(0, 0, 90, 3, φ = 90, fcol=:black)
plotpolygon2(0, 0, 90, 3, φ = 30, fcol=:black)
for θ = 0:30:330
    plotline(0, 0, 100*cos(radian(θ)), 100*sin(radian(θ)), col=:white, lwd=2)
end
newx = 200
plotcircle(newx, 0, 90, fcol=:black)
plotcircle(newx, 0, 82, fcol=:white)
plotpolygon2(newx, 0, 78, 3, φ = 90, fcol=:black)
plotpolygon2(newx, 0, 78, 3, φ = 30, fcol=:black)
for θ = 0:30:330
    plotline(newx, 0, newx+81*cos(radian(θ)), 81*sin(radian(θ)), col=:white, lwd=2)
end
plotend()

Julia に翻訳--024 連関比率法

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

連関比率法
http://aoki2.si.gunma-u.ac.jp/R/SeasonalIndex.html

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

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

label 指定は毎度毎度 label="" をやることが多いので煩わしいが，位置指定をしなくても自動的にやってくれるのがときにはうれしい。

==========#

using Statistics, Plots

function seasonalindex(x; xlab="", ylab="", main="",
                       lty1=:dash, lty2=:solid,
                       col1=:blue, col2=:black,
                       pch1=:pentagon, pch2=:circle,
                       label1="粗データ", label2="季節調整済みデータ")
    n = length(x)
    y = reshape(x./vcat(x[1], x[1:end-1]), 4, :)
    mean1 = mean(y, dims=2)
    mean2 = mean1 ./ exp(mean(log.(mean1)))
    mean2[1] = 1
    chainindex = cumprod(mean2, dims=1)
    sindex = chainindex ./ mean(chainindex)
    z = vec(reshape(x, 4, :)./sindex)
    pyplot()
    plt = plot(1:n, x, linestyle=:dash, linecolor=col1, label=label1,
               xlabel=xlab, ylabel=ylab, title=main)
    scatter!(1:n, x, markercolor=col1, markerstrokecolor=col1,
             markershape=pch1, label="")
    plot!(1:n, z, linestyle=lty2, linecolor=col2, label=label2)
    scatter!(1:n, z, markercolor=col2, markerstrokecolor=col2,
             markershape=pch2, label="")
    display(plt)
    return Dict(:seasonalindex => sindex, :correcteddata => z)
end

AirPassengers = [
    112 118 132; 129 121 135; 148 148 136; 119 104 118; 115 126 141;
    135 125 149; 170 170 158; 133 114 140; 145 150 178; 163 172 178;
    199 199 184; 162 146 166; 171 180 193; 181 183 218; 230 242 209;
    191 172 194; 196 196 236; 235 229 243; 264 272 237; 211 180 201;
    204 188 235; 227 234 264; 302 293 259; 229 203 229; 242 233 267;
    269 270 315; 364 347 312; 274 237 278; 284 277 317; 313 318 374;
    413 405 355; 306 271 306; 315 301 356; 348 355 422; 465 467 404;
    347 305 336; 340 318 362; 348 363 435; 491 505 404; 359 310 337;
    360 342 406; 396 420 472; 548 559 463; 407 362 405; 417 391 419;
    461 472 535; 622 606 508; 461 390 432];

AP = sum(AirPassengers, dims=2);
seasonalindex(AP, xlab="期", ylab="人数")
Dict{Symbol,Array{Float64,N} where N} with 2 entries:
  :correcteddata => [387.503, 376.253, 369.431, 390.518, 408.912, …
  :seasonalindex => [0.934187; 1.02325; 1.16936; 0.8732]

Julia に翻訳--023 パレート図

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

パレート図
http://aoki2.si.gunma-u.ac.jp/R/Pareto.html

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

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

縦軸を 2 つ使うのは，Julia ではできない状態が長く続いているようだ。しかたないので自分で描く。

==========#

using Plots

function pareto(x; name="", ymax=sum(x), sortflag=true, col=:blue, lwd=1,
                main="", xlab="", ylab="度数",  ylab2="累積％")
    !sortflag || sort!(x, rev=true)
    pyplot()
    n = length(x)
    total = sum(x)
    plt = bar(x, bar_width=1, ylims=(0, ymax), color=col, alpha=0.2,
              grid=false, tick_direction=:out, xlims=(0, n+2),
              xlabel=xlab, ylabel=ylab, title=main, label="")
    plot!(1.5:n+0.5, cumsum(x), linewidth=lwd, label="")
    xticks!(1:n, name)
    # 右の軸を描く
    plot!([n+1, n+1], [0, total], color=:black, label="")
    for i = 0:0.1:1
        plot!([n+1, n+1.05], [i*total, i*total], color=:black, label="")
        annotate!(n+1.05, i*total, text(" "*string(Int(100i)), 8, :left))
    end
    annotate!(n+1.8, total/2, text(ylab2, rotation=90, 11))
    display(plt)
end

x = [56, 15, 38, 8, 10, 4, 2, 2, 1, 1];
name = ["つや消え", "気泡", "異物", "ふくれ", "すりきず", "汚れ",
    "割れ", "打ちきず", "凹凸", "ながれ"];
pareto(x, name=name, xlab="樹脂部品の不良原因")