裏 RjpWiki

文字通り,RjpWiki の裏を行きます
R プログラム コンピュータ・サイエンス 統計学

Julia に翻訳--029

2021年03月05日 | ブログラミング

#==========
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 に翻訳--027

2021年03月05日 | ブログラミング

#==========
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 に翻訳--028

2021年03月05日 | ブログラミング

#==========
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 に翻訳--026

2021年03月05日 | ブログラミング

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

2021年03月05日 | ブログラミング

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

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

ファイル名: plot.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], [y1, y1, y2, y2], linecolor=col,
              linestyle=lty, linewidth=lwd, seriestype=:shape, label="")
    else
        plot!([x1, x2, x2, x1], [y1, y1, y2, y2], 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

2021年03月05日 | ブログラミング

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

2021年03月05日 | ブログラミング

#==========
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="樹脂部品の不良原因")

コメント

Julia に翻訳--022

2021年03月04日 | ブログラミング

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

ローレンツ曲線(ジニ係数)
http://aoki2.si.gunma-u.ac.jp/R/Gini-index.html

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

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

特に難しいところはなかった。

==========#

using Plots

function giniindex(y, main="", xlab="", ylab="")
    all(y .>= 0) || error("all y must be positive")
    n = length(y)
    sort!(y)
    y = cumsum(y)
    y = vcat(0, y ./ y[end])
    x = range(0, 1, length=n+1)
    pyplot()

    p = plot(x, y, color=:blue, grid=false, tick_direction=:out,
             aspect_ratio=1, size=(600, 600), label="",
             title=main, xlabel=xlab, ylabel=ylab)
    plot!([0, 1, 1, 0, 0, 1], [0, 0, 1, 0, 1, 1], color=:black, label="")
    display(p)
    return 2sum(x .- y) / n
end

x = [3, 4, 4, 5, 5, 5, 5, 6, 6, 7]
giniindex(x)

コメント

Julia に翻訳--021

2021年03月04日 | ブログラミング

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

ステム・アンド・リーフ
http://aoki2.si.gunma-u.ac.jp/R/stem.html

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

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

度数分布表を作るだけで FreqTables, printf を使うだけで Printf を読むのは大げさだなあ

==========#

using FreqTables, Printf

function stemandleaf(d, f=-99)
    function getfactor(x, minx, maxx)
        for i in vcat(1:10, -10:-1)
            ll = floor(Int, maxx*10.0^i)-floor(Int, minx*10.0^i)
            if 2 <= ll < 19
                 return 10.0^i
            end
        end
        return 1
    end

    DUMMY = 99
    MINUS = -0.1
    f = f == -99 ? getfactor(d, minimum(d), maximum(d)) : 10^f
    temp = floor.(Int, d*f*10)
    stem = floor.(Int, temp/10)
    leaf = abs.(temp) - abs.(10stem)
    stem = stem == 0 ? (d > 0 ? 0 : MINUS) : stem
    minstem = minimum(stem)
    maxstem = maximum(stem)
    stem = vcat(stem, minstem:maxstem)
    leaf = vcat(leaf, fill(DUMMY, maxstem-minstem+1))
    if maxstem > 0 && minstem < 0
        stem = vcat(stem, MINUS)
        leaf = vcat(leaf, DUMMY)
    end
    res = freqtable(stem, leaf)
    rownames = names(res)[1]
    colnames = names(res)[2]
    nr, nc = size(res)
    for i = 1:nr
        stemtemp = rownames[i]
        @printf("%5d|", stemtemp == MINUS ? "-0" : stemtemp)
        for le = 1:nc
            if names(res)[2][le] != DUMMY
                print(string(colnames[le]) ^ res[i, le])
            end
        end
        print("\n")
    end
    println("stem * $(1/f)")
end

using Random
Random.seed!(123);
d = floor.(Int, randn(200) .* 20 .+ 100);
stemandleaf(d)

    3|7
    4|11
    5|4
    6|14457999
    7|01123346777888899
    8|0011122223444455556667788999
    9|00011122233333344445555566666677777899
   10|000000011222334444444555566666677778888999
   11|0011123333445566666778899
   12|00011223335556689
   13|00012333468
   14|00125566
   15|25
stem * 10.0

コメント

Julia に翻訳--020

2021年03月04日 | ブログラミング

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

シンプレックス法によるパラメータ推定
http://aoki2.si.gunma-u.ac.jp/R/simplex.html

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

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

loops は for ループ中のローカル変数なので,ループが回りきったかどうかを知るために for ループの外で参照しようとしてもできない(loops は存在しないというエラーメッセージが出る)。
その他にも,元の R プログラムでは求めるパラメータもローカル変数なので for ループの外で利用することができないなどのため,パラメータが求まったときの処理を for ループ内で行うように変更した。

==========#

using Plots

function simplex(fun, start, x, y; MAXIT = 10000, EPSILON = 1e-7,
                 LO = 0.8, HI = 1.2, plotflag = false)
    # one line function definition
    residual(x, y, p) = sum((y .- fun(x, p)) .^ 2)

    ip = length(start)
    ip1 = ip + 1
    ip2 = ip + 2
    ip3 = ip + 3
    pa = reshape(repeat(start, ip3), ip, :)
    for i = 1:ip
        pa[i, i] = start[i] * rand(1)[1] * (HI - LO) + LO
    end
    res = vcat([residual(x, y, pa[:, i]) for i = 1:ip1], 0, 0)
    converge = false
    for loops = 1:MAXIT
        res0 = res[1:ip1]
        mx = argmax(res0)
        mi = argmin(res0)
        s = sum(pa[:, 1:ip1], dims = 2)
        if res[mx] < EPSILON || res[mi] < EPSILON ||
           (res[mx] - res[mi]) / res[mi] < EPSILON
            converge = true
            parameters = pa[:, mi]
            residuals = residual(x, y, parameters)
            if plotflag
                pyplot()
                plt = scatter(x, y, tick_direction = :out, grid = false,
                    markercolor = :blue, label = "")
                x = range(minimum(x), maximum(x), length = 1000)
                plot!(x, fun(x, parameters), label = "")
            end
            display(plt)
            return Dict(:converge => converge, :parameters => parameters, :residuals => residuals)
        end
        i = ip2
        pa[:, ip2] = (2 * s - ip2 * pa[:, mx]) / ip
        res[ip2] = residual(x, y, pa[:, ip2])
        if res[ip2] < res[mi]
            pa[:, ip3] = (3 * s - (2 * ip1 + 1) * pa[:, mx]) / ip
            res[ip3] = residual(x, y, pa[:, ip3])
            if res[ip3] <= res[ip2]
                i = ip3
            end
        elseif res[ip2] > res[mx]
            pa[:, ip3] = s / ip1
            res[ip3] = residual(x, y, pa[:, ip3])
            if res[ip3] >= res[mx]
                for i = 1:ip1
                    if i != mi
                        pa[:, i] = (pa[:, i] + pa[:, mi]) * 0.5
                        res[i] = residual(x, y, pa[:, i])
                    end
                end
                i = 0
            else
                i = ip3
            end
        end
        if i > 0
            pa[:, mx] = pa[:, i]
            res[mx] = res[i]
        end
    end
    println("not converged!")
end

x = 1:10; # x 値
y = [3, 8, 15, 35, 57, 80, 92, 95, 99, 100]; # y 値

# あてはめるモデル関数
fun(x, p) = p[1] ./ (1 .+ p[2] .* exp.(-p[3] .* x))

simplex(fun, [80, 70, 0.5], x, y, plotflag = true)

コメント

Julia に翻訳--018

2021年03月03日 | ブログラミング

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

S の plot.design 関数
http://aoki2.si.gunma-u.ac.jp/R/plot-design.html

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

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

! で重ね描きできるのは便利。

==========#

using DataFrames, Statistics, Plots

function plotdesign(data; FUN=mean)
    nv = size(data, 2)
    pyplot()
    name = names(data)
    p = plot(xlims=(1.5, nv+0.5), ylabel=name[1], yguidefontsize=14,
             tick_direction=:out, grid=false, label="")
    minx, maxx = Inf, -Inf
    for i in 2:nv
        gd = groupby(data, i)
        a = combine(gd, 1 => FUN)
        x = repeat([i], size(a, 1))
        plot!(x, a[:, 2], color=:black, label="")
        scatter!(x, a[:, 2], color=:red, label="")
        for j = 1:length(x)
            minx = min(minx, a[j, 2])
            maxx = max(maxx, a[j, 2])
            annotate!(x[j], a[j, 2], text(" " * a[j, 1], 12, :left), label="")
        end
    end
    w = (maxx - minx) * 25/400
    xticks!(2:nv, name[2:nv], xtickfontsize=14,
            ylims=(minx - w, maxx + w))
end

using Random;
Random.seed!(123);
n = 50
value = randn(n) .* 10 .+ 100;
gender = rand(["male", "female"], n);
area = rand(["rural", "urban"], n);
density = rand(["hi", "med", "lo"], n);
dist = rand(["a", "b", "c", "d", "e"], n);
# 第1列に連続変数,2列目以降にカテゴリー変数
data = DataFrame(:value => value, :gender => gender, :area => area,
                 :density => density, :dist => dist);
plotdesign(data)


plotdesign(data[:, 1:3])

コメント

Julia に翻訳--019

2021年03月03日 | ブログラミング

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

散布図,確率楕円,回帰直線,信頼限界帯,MA regression,RMA regression
http://aoki2.si.gunma-u.ac.jp/R/scatter.html

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

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

あれこれオプションを追加すると,引数も増える。

==========#

using Statistics, Distributions, LinearAlgebra, Plots

function scatterplot(x, y;
                     el =false, el2 =(:black,   :solid, 0.5),
                     lrl=false, lrl2=(:black,   :solid, 0.5),
                     cb =false, cb2 =(:blue,    :dot,   0.5),
                                cb3 =(:brown,   :dash,  0.5),
                     ma =false, ma2 =(:magenta, :solid, 0.5),
                     rma=false, rma2=(:gray,    :solid, 0.5),
                     marker=(:red, :red, 0.25), factor=0.15,
                     α=0.05, acc=2000, xlab="", ylab="")
    function abline(a, b, col, sty, wid)
        plot!([minx, maxx], [a + b * minx, a + b * maxx],
              linecolor=col, linestyle=sty, linewidth=wid, label="")
    end

    function ellipsedraw(draw, el2, α, acc)
        λ = eigvals(v)
        a = sqrt(vxy^2 / ((λ[1] - vx)^2 + vxy^2))
        b = (λ[1] - vx) * a / vxy
        θ = atan(a/b)
        k = sqrt(-2log(α))
        l1 = sqrt(λ[2]) * k
        l2 = sqrt(λ[1]) * k
        x2 = range(-l1, l1, length=acc)
        tmp = [tmp0 < 0 ? 0 : tmp0 for tmp0 in 1 .- x2 .^ 2 ./ l1 ^ 2]
        y2 = l2 .* sqrt.(tmp)
        x2 = vcat(x2, reverse(x2))
        y2 = vcat(y2, -reverse(y2))
        s0 = sin(θ)
        c0 = cos(θ)
        xx = c0 .* x2 .+ s0 .* y2 .+ mean(x)
        yy = -s0 .* x2 .+ c0 .* y2 .+ mean(y)
        if draw
            plot!(xx, yy, linecolor=el2[1], linestyle=el2[2], linewidth=el2[3], label="")
        end
        return minimum(xx), maximum(xx), minimum(yy), maximum(yy)
    end

    function conflimit(lrl2, cb, cb2, cb3, α)
        b = vxy / vx
        a = mean(y) - b * mean(x)
        abline(a, b, lrl2[1], lrl2[2], lrl2[3])
        if cb
            sx2 = vx*(n-1)
            x1 = range(minx, maxx, length=2000)
            y1 = a .+ b .* x1
            ta = cquantile(TDist(n - 2), α/2)
            Ve = (vy - vxy^2 / vx) * (n - 1) / (n - 2)
            temp = ta * sqrt(Ve) * sqrt.(1/n .+ (x1 .- mean(x)) .^ 2 ./ sx2)
            y2 = y1 - temp;
            plot!(x1, y2, linecolor=cb2[1], linestyle=cb2[2], linewidth=cb2[3], label="")
            y2 = y1 + temp;
            plot!(x1, y2, linecolor=cb2[1], linestyle=cb2[2], linewidth=cb2[3], label="")
            temp = ta * sqrt(Ve) * sqrt.(1 + 1/n .+ (x1 .- mean(x)) .^ 2 ./ sx2);
            y2 = y1 - temp;
            plot!(x1, y2, linecolor=cb3[1], linestyle=cb3[2], linewidth=cb3[3], label="")
            y2 = y1 + temp;
            plot!(x1, y2, linecolor=cb3[1], linestyle=cb3[2], linewidth=cb3[3], label="")
        end
        print("LS(least squares) ")
        println("intercept = $(round(a, digits=5)), slope = $(round(b, digits=5))")
    end

    function MA(ma2)
        b = vxy / (eigvals(v)[2] - vy)
        a = mean(y) - b * mean(x)
        abline(a, b, ma2[1], ma2[2], ma2[3])
        print("MA(major axis) ")
        println("intercept = $(round(a, digits=5)), slope = $(round(b, digits=5))")
        Dict(:intercept => a, :slope => b)
    end

    function RMA(rma2)
        b = sign(vxy) * sqrt(vy / vx)
        a = mean(y) - b * mean(x)
        abline(a, b, rma2[1], rma2[2], rma2[3])
        print("RMA(reduced major axis) ")
        println("intercept = $(round(a, digits=5)), slope = $(round(b, digits=5))")
        Dict(:intercept => a, :slope => b)
    end

    function defminmax(minx0, maxx0, x)
        margin = (max(maxx0, maximum(x)) - min(minx0, minimum(x))) * factor
        return minx0 - margin, maxx0 + margin
    end

    pyplot()
    v = cov(hcat(x, y))
    vx = v[1, 1]
    vy = v[2, 2]
    vxy = v[1, 2]
    minx, maxx, miny, maxy = ellipsedraw(false, el2, α, acc)
    minx, maxx = defminmax(minx, maxx, x)
    miny, maxy = defminmax(miny, maxy, y)
    n = length(x)
    p = scatter(x, y, xlims=(minx, maxx), ylims=(miny, maxy),
                grid=false, tick_direction=:out, color=marker[1],
                markerstrokecolor=marker[2], alpha=marker[3],
                xlabel=xlab, ylabel=ylab, label="")
    !el  || ellipsedraw(true, el2, α, acc)
    !lrl || conflimit(lrl2, cb, cb2, cb3, α)
    !ma  || MA(ma2)
    !rma || RMA(rma2)
    display(p)
end

x = [132, 146, 140, 196, 132, 154, 154, 168, 140, 140, 156, 114, 134, 116, 150, 178, 150, 120, 150, 146];
y = [90, 90, 84, 96, 90, 90, 74, 92, 60, 82, 80, 62, 80, 80, 76, 98, 86, 70, 80, 80];
scatterplot(x, y, el=true, lrl=true, cb=true)


z = randn((1000, 20));
x2 = mean(z[:, 1:12], dims=2);
y2 = mean(z[:, 8:20], dims=2);
scatterplot(x2, y2, el=true)

コメント

Julia に翻訳--017

2021年03月03日 | ブログラミング

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

ROC 曲線と ROC 曲線下面積
http://aoki2.si.gunma-u.ac.jp/R/ROC.html

ファイル名: roc.jl  関数名: roc (同じ名前で 2 個の関数)

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

引数の個数が違うので,同じ名前で定義しても混乱がない。
度数分布の階級は,R の pretty() を使う。

==========#

using Plots, DataFrames, Printf, RCall

# for raw data
function roc(disease, normal; lowest=Inf, width=Inf, verbose=true)
    function myhist(x, brks)
        k = length(brks)
        freq = zeros(Int, k)
        for i in 1:k-1
            freq[i] = sum(brks[i] .<= x .< brks[i+1])
        end
        freq[k] = sum(x .>= brks[k])
        return freq
    end
    x = vcat(disease, normal);
    minimum(x)
    R"prettynumber = pretty($x)"
    @rget prettynumber
    minx = minimum(x)
    maxx = maximum(x)
    if verbose
        println("最小値 = $minx")
        println("最大値 = $maxx")
    end
    if lowest == Inf || width == Inf
        R"prettynumber = pretty($x, 10)"
        @rget prettynumber
        lowest = prettynumber[1]
        width = diff(prettynumber)[1]
        if verbose
            println("最小値より小さいキリのよい数値 = $lowest")
            println("度数分布を作成する階級幅の切りのよい数値 = $width")
        end
    end
    brks = collect(range(lowest, maxx+width, step=width))
    roc(float.(brks), myhist(disease, brks), myhist(normal, brks), verbose=verbose)
end

# for frequency data
function roc(x, disease, normal; verbose=true)
    k = length(x)
    (k == length(disease) == length(normal)) || error("length is not same")
    sensitivity = vcat(reverse(cumsum(reverse(disease)))/sum(disease), 0)
    falsepositiverate = vcat(reverse(cumsum(reverse(normal)))/sum(normal), 0)
    specificity = 1 .- falsepositiverate
    pyplot()
    p = plot(falsepositiverate, sensitivity, color=:black, grid=false,
             tick_direction=:out, aspect_ratio=1, label="")
    scatter!(falsepositiverate, sensitivity, color=:black, label="")
    plot!([0, 1, 1, 0, 0], [0, 0, 1, 1, 0], color=:black,
          linewidth=0.5, label="")
    display(p)
    cindex = sum((falsepositiverate[i]-falsepositiverate[i+1]) *
                 (sensitivity[i+1]+sensitivity[i])/2 for i = 1:k)
    result = DataFrame("Value" => x, "Disease" => disease, "Normal" => normal,
                        "Sensitivity" => sensitivity[1:end-1],
                        "Specificity" => specificity[1:end-1],
                        "F.P. rate" => falsepositiverate[1:end-1])
    if verbose
        @printf("%3s  %11s  %7s  %6s  %11s  %11s  %9s\n",
                "", "Value", "Disease", "Normal", "Sensitivity",
                "Specificity", "F.P. rate")
        for i = 1:k
            @printf("%3d  %11.5g  %7d  %6d  %11.3f  %11.3f  %9.3f\n",
                    i, x[i], disease[i], normal[i], sensitivity[i],
                    specificity[i], falsepositiverate[i])
        end
        @printf("c index = %.5g\n", cindex)
    end
    return Dict("result" => result, "cindex" => cindex)
end

FRA = [100, 220, 230, 240, 250, 260, 270, 280, 290, 300, 320, 340, 360, 400]
FRAdisease = [3, 2, 1, 4, 7, 4, 16, 5, 3, 9, 10, 5, 10, 21]
FRAnormal  = [25, 7, 19, 17, 7, 8, 7, 6, 2, 2, 0, 0, 0, 0]
roc(FRA, FRAdisease, FRAnormal)

           Value  Disease  Normal  Sensitivity  Specificity  F.P. rate
  1          100        3      25        1.000        0.000      1.000
  2          220        2       7        0.970        0.250      0.750
  3          230        1      19        0.950        0.320      0.680
  4          240        4      17        0.940        0.510      0.490
  5          250        7       7        0.900        0.680      0.320
  6          260        4       8        0.830        0.750      0.250
  7          270       16       7        0.790        0.830      0.170
  8          280        5       6        0.630        0.900      0.100
  9          290        3       2        0.580        0.960      0.040
 10          300        9       2        0.550        0.980      0.020
 11          320       10       0        0.460        1.000      0.000
 12          340        5       0        0.360        1.000      0.000
 13          360       10       0        0.310        1.000      0.000
 14          400       21       0        0.210        1.000      0.000
c index = 0.88215

using Random
Random.seed!(123);

disease = floor.(Int, randn(200) .* 10 .+ 110);
normal = floor.(Int, randn(300) .* 10 .+ 100);
roc(disease, normal)

最小値 = 72
最大値 = 137
最小値より小さいキリのよい数値 = 70
度数分布を作成する階級幅の切りのよい数値 = 5
           Value  Disease  Normal  Sensitivity  Specificity  F.P. rate
  1           70        0       5        1.000        0.000      1.000
  2           75        1       3        1.000        0.017      0.983
  3           80        2       9        0.995        0.027      0.973
  4           85        1      33        0.985        0.057      0.943
  5           90        8      48        0.980        0.167      0.833
  6           95       17      63        0.940        0.327      0.673
  7          100       28      46        0.855        0.537      0.463
  8          105       38      37        0.715        0.690      0.310
  9          110       42      25        0.525        0.813      0.187
 10          115       25      19        0.315        0.897      0.103
 11          120       17       8        0.190        0.960      0.040
 12          125       11       4        0.105        0.987      0.013
 13          130        8       0        0.050        1.000      0.000
 14          135        2       0        0.010        1.000      0.000
 15          140        0       0        0.000        1.000      0.000
c index = 0.75928

コメント

Julia に翻訳--016

2021年03月02日 | ブログラミング

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

チャーノフの顔グラフ
http://aoki2.si.gunma-u.ac.jp/R/face.html

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

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

変数スコープに嵌まってしまった。

==========#

using Plots

function face(x, size0 = 480)

    function arc1(x1, y1, r, l)
        sign2 = l > 0 ? -1 : 1
        θ = sign2 * acos(x1 / r)
        y1 -= sign2 * sqrt(r^2 - x1^2)
        if l <= 0
            arc(0, y1, r, θ, π - θ)
        else
            arc(0, y1, r, π - θ, 2π + θ)
        end
    end

    function arc(ox, oy, r, θstart, θend)
        step2 = min(0.1, (θend - θstart) * 0.1)
        interval = θstart:step2:θend
        plot!(r .* cos.(interval) .+ ox, r .* sin.(interval) .+ oy,
              color = :black, label = "")
    end

    function ellipse(ox, oy, ra, rb, θaxis, θstart, θend)
        θend += 2π * (θend <= θstart)
        temp1 = ra * rb
        temp2 = 30 / (ra + rb)
        k = floor(Int, (θend - θstart) / temp2) + 2
        xx = zeros(k)
        yy = zeros(k)
        for i = 1:k-1
            factor = temp1 / sqrt((ra * sin(θstart))^2 + (rb * cos(θstart))^2)
            xx[i] = factor * cos(θaxis + θstart)
            yy[i] = factor * sin(θaxis + θstart)
            θstart += temp2
        end
        factor = temp1 / sqrt((ra * sin(θend))^2 + (rb * cos(θend))^2)
        xx[k] = factor * cos(θaxis + θend)
        yy[k] = factor * sin(θaxis + θend)
        plot!(ox .+ xx, oy .+ yy, color = :black, label = "")
    end

    pi2 = 2π
    p = plot(xlims = (-500, 500), ylims = (-500, 500), lab = "", ticks = false,
            showaxis = false, aspect_ratio = 1, label = "" )
    size2 = size0 * (1 + x[1]) / 2
    θ = (π / 4) * (2 * x[2] - 1)
    h = size0 * (1 + x[3]) / 2
    x1 = size2 * cos(θ)
    y1 = size2 * sin(θ)
    ak = 1 - x[4]^2
    oy1 = (ak * x1^2 + y1^2 - h^2) / (2 * (y1 - h))
    rb1 = h - oy1
    ra1 = rb1 / sqrt(ak)
    θstart = atan((y1 - oy1) / x1)
    θend = π - θstart
    ellipse(0, oy1, ra1, rb1, 0, θstart, θend)
    ak = 1 - x[5]^2
    oy2 = (ak * x1^2 + y1^2 - h^2) / (2 * (y1 + h))
    rb2 = h + oy2
    ra2 = rb2 / sqrt(ak)
    θend = atan((y1 - oy2) / x1)
    θstart = π - θend
    ellipse(0, oy2, ra2, rb2, 0, θstart, θend)
    y = h * x[6]
    plot!([0, 0], [y, -y], color = :black, label = "")
    pm = -h * (x[7] + (1 - x[7]) * x[6])
    wm = sqrt(ra2^2 * (1 - (pm - oy2)^2 / rb2^2))
    if x[8] == 0
        plot!([-wm / 2, wm / 2], [pm, pm], color = :black, label = "")
    else
        r = h / abs(x[8])
        am = x[9] * r
        x1 = ifelse(am > wm, x[9] * wm, am)
        l = ifelse(x[8] <= 0, -1, 1)
        y1 = pm - l * (r - sqrt(r^2 - x1^2))
        arc1(x1, y1, r, l)
    end
    ye = h * (x[10] + (1 - x[10]) * x[6])
    we = sqrt(ra1^2 * (1 - (ye - oy1)^2 / rb1^2))
    xe = we * (1 + 2 * x[11]) / 4
    θ = (2 * x[12] - 1) * π / 5
    ra3 = x[14] * min(xe, we - xe)
    rb3 = sqrt(ra3^2 * (1 - x[13]^2))
    ellipse(xe, ye, ra3, rb3, θ, 0, pi2)
    ellipse(-xe, ye, ra3, rb3, π - θ, 0, pi2)
    re = ra3 / sqrt(cos(θ)^2 + sin(θ)^2 / x[13]^2)
    shift = re * (2 * x[15] - 1)
    arc(xe - shift, ye, 3, 0, pi2)
    arc(-xe - shift, ye, 3, 0, pi2)
    θ2 = 2 * (1 - x[17]) * (π / 5)
    θ3 = θ >= 0 ? θ + θ2 : θ - θ2
    len = re * (2 * x[18] + 1) / 2
    x0 = len * cos(θ3)
    x1 = xe .- [x0, -x0]
    y0 = len * sin(θ3)
    y1 = ye + 2 * (x[16] + 0.3) * ra3 * x[13] .- [y0, -y0]
    plot!(x1 .- shift, y1, color = :black, label = "")
    plot!(-x1 .- shift, y1, color = :black, label = "")
    #display(p)
    return p
end

function facedata(d; pos = collect(1:18))
    lo = [0.2, 0.2, 0.2, 0.2, 0.2, 0.1, 0.2, -5, 0.2,
          0.1, 0.1, 0.3, 0.1, 0.3, 0.1, 0.1, 0.1, 0.1]
    hi = [0.8, 0.8, 1, 0.8, 0.8, 0.4, 0.8, 5, 0.8, 0.7,
          0.9, 0.7, 0.9, 0.9, 0.9, 0.9, 1, 0.9]
    fx = [0.5, 0.5, 1, 0.5, 0.5, 0.2, 0.5, 0, 0.5, 0.4,
          0.5, 0.5, 0.5, 0.6, 0.5, 0.5, 1, 0.5]
    n = size(d, 1)
    mind = vec(minimum(d, dims = 1))
    maxd = vec(maximum(d, dims = 1))
    local x = zeros((n, 18))
    for i = 1:18 # i = 3
        k = pos[i]
        if k > 0
            mindk = mind[k]
            maxdk = maxd[k]
            if mindk == maxdk
                x[:, i] = fx[i]
            else
                x[:, i] .= (d[:, k] .- mindk) ./ (maxdk - mindk) .* (hi[i] - lo[i]) .+ lo[i]
            end
        end
    end
    return x
end

nv = 18
using Random
Random.seed!(123);
d = randn((9, nv));
a = facedata(d);
pyplot()
a1 = face(a[1, :])
a2 = face(a[2, :])
a3 = face(a[3, :])
a4 = face(a[4, :])
plot(a1, a2, a3, a4, layout = (2, 2))

コメント

Julia に翻訳--015

2021年03月02日 | ブログラミング

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

レーダーチャート
http://aoki2.si.gunma-u.ac.jp/R/radar.html

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

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

条件 || 条件を満たさないときに実行される式
条件 && 条件を満たすときに実行される式

は,Julia 独特だが,ちょっとわかりにくい。

==========#

using Statistics, Plots

function radar(df; maxx=Inf, minx=Inf, normalize=true, col="", lty="", title="")
    function drawnet(x, border, lty)
        scale = (x .- minx) ./ (maxx .- minx)
        scale = vcat(scale, scale[1])
        if all(0 .< scale .<= 1.2)
            plot!(vcat(cosine, cosine[1]).*scale, vcat(sine, sine[1]).*scale,
                  seriestype=:path, label="", linecolor=border, linestyle=lty)
        end
    end
    pyplot()
    axisname = names(df)
    df = Matrix(df)
    n, m = size(df)
    m >= 3 || error("3変数以上であるべし!")
    p = plot(xlimit=(-1.4, 1.4), ylimit=(-1.2, 1.2), aspect_ratio=1, axis=false,
             ticks=false, xlabel="", ylabel="", title=title, size=(600, 400))
    θ = π .*(0.5 .- collect(0:m-1) .* 2 ./ m)
    cosine = cos.(θ)
    sine = sin.(θ)
    if normalize
        means = mean(df, dims=1)
        stds = std(df, dims=1)
        for i = 1:m
            df[:, i] = (df[:, i] .- means[i]) ./ stds[i]
        end
        maxx = maximum(df)
        minx = minimum(df)
        w = (maxx - minx) * 0.1
        maxx = maxx + w
        minx = minx - w
        for i = -3:3
            x = fill(i, m)
            drawnet(x, :gray, i==0 ? :solid : :dot)
        end
    else
        maxx != Inf || (maxx = vec(maximum(df, dims=1)))
        minx != Inf || (minx = vec(minimum(df, dims=1)))
        w = (maxx - minx) * 0.1
        maxx = maxx + w
        minx = minx - w
        for i = 1:5
            x = cosine * i / 5
            y = sine * i / 5
            plot!(vcat(x, x[1]), vcat(y, y[1]), seriestype=:path,
                    linecolor=:gray, linestyle=:dot, label="")
        end
    end
    for i = 1:m
        plot!([0, cosine[i]*1.05], [0, sine[i]*1.05], seriestype=:path,
                color=:gray, label="")
        annotate!(cosine[i]*1.1, sine[i]*1.1, text(axisname[i], 10,
                i > m/2 ? :right : :left))
    end
    col != "" || (col = fill(:red, m))
    lty != "" || (lty = fill(:solid, m))
    for i = 1:n
        drawnet(df[i, :], col[i], lty[i])
    end
    display(p)
end

using RDatasets
swiss = dataset("datasets", "swiss")
test = swiss[1:5, 2:7]
radar(test)


radar(test, col=[:black, :red, :green4, :blue, :brown])


radar(test, maxx=[100, 50, 20, 20, 100, 25], normalize=false,
        col=[:black, :red, :green4, :blue, :brown])


radar(test, maxx=[100, 50, 20, 20, 100, 25],
        minx=[70, 10, 0, 0, 0, 20], normalize=false,
        col=[:black, :red, :green4, :blue, :brown])


radar(test, maxx=100, minx=0, normalize=false,
        col=[:black, :red, :green4, :blue, :brown])

コメント