裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

Julia に翻訳--019 散布図,確率楕円,回帰直線,信頼限界帯,MA regression,RMA regression

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)

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia に翻訳--018 S の plot.design 関数

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])

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia に翻訳--017 ROC 曲線,ROC 曲線下面積

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

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村