裏 RjpWiki

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

Julia に翻訳--147　総当たり法による重回帰分析

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

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

総当たり法による重回帰分析
http://aoki2.si.gunma-u.ac.jp/R/All_possible_subset_selection.html

ファイル名: allpossiblesubsetselection.jl
関数名:    allpossiblesubsetselection(dat; vname=[], limit = 10)
          printallpossiblesubsetselection(obj; models = 20, sortby = "adj")

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

分析結果を並べ替えて表示する関数を別立てにしている。

==========#


using Statistics, Combinatorics, Printf

function allpossiblesubsetselection(dat; vname=[], limit = 10)
    nv = size(dat, 2)
    nv <= limit || error("独立変数が $limit 個以上である（多すぎる）。\nlimit 引数で変更できる。")
    length(vname) != 0 || (vname = vcat(["x" * string(i) for i = 1:nv-1], "y"))
    nv -= 1
    n = 2 ^ nv
    depname = vname[end]
    indepname = vname[1:end-1]
    indepdat = copy(dat[:, 1:nv])
    formula = fill("", n);
    results = zeros(n, 4);
    for (i, j) in enumerate(Iterators.product([[0, 1] for i = 1:nv]...))
        suf = [k == 1 for k in j]
        str = indepname[suf]
        dat1 = hcat(indepdat[:, suf], dat[:, end])
        formula[i] = depname * " ~ " * join(str, " + ")
        results[i, :] .= mreg(dat1)
    end
    res = Dict(:results => (results[2:end, :], formula[2:end]))
    printallpossiblesubsetselection(res; models = 20, sortby = "rsq")
    res
end

function mreg(dat) # 回帰分析を行い R2, R2a, AIC, Loglik を返す
    n, nc = size(dat)
    nv = nc - 1
    r = cor(dat)
    betas = r[1:end-1, 1:end-1] \ r[1:end-1, end]
    Sr = sum(r[1:nv, nc] .* betas)
    Se = 1 - Sr
    Vy = var(dat[:, nc]) * (n - 1)
    R2 = Sr
    R2a = 1 - Se / (n - nv - 1) * (n - 1)
    LogLik = -0.5n * (log(2π) + 1 - log(n) + log(Se * Vy))
    AIC = 2nc + 2 - 2LogLik
    [R2, R2a, AIC, LogLik]
end

function printallpossiblesubsetselection(obj; models = 20, sortby = "adj") # "rsq", "AIC", "loglik"
    formula = obj[:results][2]
    results = obj[:results][1]
    key = indexin([sortby],["adj", "rsq", "AIC", "loglik"])[1]
    o = sortperm(results[:, key], rev=(sortby!="AIC"))
    results = results[o, :]
    formula = formula[o]
    @printf("\n R square  Adjusted R square       AIC      log-likelihood     Formula\n")
    for i in 1:min(models, size(results, 1))
        @printf(" %8.5f %13.5f     %13.5f  %13.5f       %s\n",
        results[i, 1], results[i, 2], results[i, 3], results[i, 4], formula[i])
    end
end

using RDatasets
iris = dataset("datasets", "iris");
vname=["SepalLength", "SepalWidth", "PetalLength", "PetalWidth"];
dat = Array(iris[:, 1:4]);
res = allpossiblesubsetselection(dat, vname=vname)
#=
 R square  Adjusted R square       AIC      log-likelihood     Formula
  0.93785       0.93657         -63.50218       36.75109       PetalWidth ~ SepalLength + SepalWidth + PetalLength
  0.92975       0.92879         -47.11940       27.55970       PetalWidth ~ SepalWidth + PetalLength
  0.92902       0.92806         -45.58470       26.79235       PetalWidth ~ SepalLength + PetalLength
  0.92711       0.92662         -43.59109       24.79555       PetalWidth ~ PetalLength
  0.74293       0.73943         147.46929      -69.73464       PetalWidth ~ SepalLength + SepalWidth
  0.66903       0.66679         183.37107      -88.68553       PetalWidth ~ SepalLength
  0.13405       0.12820         327.64025     -160.82012       PetalWidth ~ SepalWidth
Dict{Symbol,Tuple{Array{Float64,2},Array{String,1}}} with 1 entry:
  :results => ([0.669028 0.666791 183.371 -88.6855; 0.134048 0.128197 327.64 -1…
=#

printallpossiblesubsetselection(res, sortby="AIC")
#=
 R square  Adjusted R square       AIC      log-likelihood     Formula
  0.93785       0.93657         -63.50218       36.75109       PetalWidth ~ SepalLength + SepalWidth + PetalLength
  0.92975       0.92879         -47.11940       27.55970       PetalWidth ~ SepalWidth + PetalLength
  0.92902       0.92806         -45.58470       26.79235       PetalWidth ~ SepalLength + PetalLength
  0.92711       0.92662         -43.59109       24.79555       PetalWidth ~ PetalLength
  0.74293       0.73943         147.46929      -69.73464       PetalWidth ~ SepalLength + SepalWidth
  0.66903       0.66679         183.37107      -88.68553       PetalWidth ~ SepalLength
  0.13405       0.12820         327.64025     -160.82012       PetalWidth ~ SepalWidth
=#

printallpossiblesubsetselection(res, sortby="loglik")
#=
 R square  Adjusted R square       AIC      log-likelihood     Formula
  0.93785       0.93657         -63.50218       36.75109       PetalWidth ~ SepalLength + SepalWidth + PetalLength
  0.92975       0.92879         -47.11940       27.55970       PetalWidth ~ SepalWidth + PetalLength
  0.92902       0.92806         -45.58470       26.79235       PetalWidth ~ SepalLength + PetalLength
  0.92711       0.92662         -43.59109       24.79555       PetalWidth ~ PetalLength
  0.74293       0.73943         147.46929      -69.73464       PetalWidth ~ SepalLength + SepalWidth
  0.66903       0.66679         183.37107      -88.68553       PetalWidth ~ SepalLength
  0.13405       0.12820         327.64025     -160.82012       PetalWidth ~ SepalWidth
=#

 

頑張ってます。クリックお願いします。

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

Julia に翻訳--146　重回帰分析，ステップワイズ変数選択

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

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

重回帰分析（ステップワイズ変数選択）
http://aoki2.si.gunma-u.ac.jp/R/sreg.html

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

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

翻訳のしがいのあるプログラムだった。

==========#

using DataFrames, Rmath, Statistics, Printf, Plots

function sreg(data; vname=[], stepwise = true, Pin = 0.05, Pout = 0.05,
              predict = false, verbose = false)

    formatpval(p) = p >= 0.0001 ? @sprintf("%.5f", p) : "< 0.0001"

    function stepout(isw)
        step += 1
        syy = ss[q1]
        lxi = lx[1:ni]
        b1 = r[q1, lxi];
        b = b1 .* sd[q1] ./ sd[lxi];
        b0 = means[q1] - sum(b .* means[lxi])
        sr = sum(b .* ss[lxi])
        se = syy - sr
        if se < 0 && abs(se / syy) < 1e-12
            se = 0
        end
        idf1 = ni
        idf2 = nc - ni - 1
        vr = sr / idf1
        ve = se / idf2
        if ve == 0
            f = fp9 = NaN
            stde = tv = pp = rep(NaN, ni)
            seb0 = tv0 = pp0 = NaN
            error("従属変数の予測値と観測値が完全に一致します。分析の指定に間違いはないですか？")
        else
            f = vr / ve
            fp9 = pf(f, idf1, idf2, false)
        end
        rhoy = se / syy
        if ve != 0
            rdiag = [r[lxi0, lxi0] for lxi0 in lxi];
            seb = rdiag .* rhoy ./ idf2 .* var[q1] ./ var[lxi];
            stde = sqrt.(seb);
            tv = abs.(b ./ stde);
            pp = pt.(tv, idf2, false) * 2;
            temp = means[lxi] ./ sd[lxi];
            seb0 = sqrt((1 / nc + transpose(temp) * r[lxi, lxi] * temp / (nc - 1)) * ve)
            tv0 = abs(b0 / seb0)
            pp0 = pt(tv0, idf2, false) * 2
        end
        append!(stde, seb0);
        append!(tv, tv0);
        append!(pp, pp0);
        adjR2 = 1 - ve / syy * (nc - 1)
        if rhoy != 0 && isw != 0
            r2change = oldrhoy - rhoy
            fchange = r2change * idf2 / rhoy
            pchange = fchange < 0 ? NaN : pf(fchange, 1, idf2, false)
        else
            pchange = NaN
        end
        if isw != 0 && verbose
            print("\n ***** ステップ $step *****   ")
            println("　　$(["編入", "除去"][isw])変数: $(vname[ip])")
        end
        multico = [r[lxi0, lxi0] for lxi0 in lxi]
        R2 = 1 - rhoy
        R = sqrt(R2)
        loglik = 0.5 * (sum( - nc * (log(2 * pi) + 1 - log(nc) + log(se))))
        AIC = 2 * ni + 4 - 2 * loglik
        results = merge(results, Dict(:lxi => lxi,
                    :variables => vname[lxi],
                    :coefficient => vcat(b, b0),
                    :standarderror => stde,
                    :tvalue => tv,
                    :tPvalue => pp,
                    :standardizedcoefficient => b1,
                    :tolerance => 1 ./ multico,
                    :VIF => multico,
                    :SS => vcat(sr, se, syy),
                    :df => vcat(ni, idf2, nc - 1),
                    :MS => vcat(vr, ve),
                    :Fvalue => f,
                    :FPvalue => fp9,
                    :R2 => R2,
                    :adjR2 => adjR2,
                    :R => R,
                    :loglik => loglik,
                    :AIC => AIC,
                    :ve => ve))
        if stepwise && !isnan(pchange)
            results = merge(results, Dict(:r2change => r2change,
                    :fchange => fchange,
                    :df1 => 1,
                    :df2 => idf2,
                    :pchange => pchange))
        end
        oldrhoy = rhoy
    end

    function fmax()
        kouho = 1:q
        if ni > 0
            suf = trues(q)
            [suf[i] = false for i in lx[1:ni]]
            kouho = (1:q)[suf]
        end
        rdiag = [r[i, i] for i in kouho]
        temp = 1 ./ (rdiag .* r[q1, q1] ./ (r[q1, kouho] .^ 2) .- 1)
        ip = argmax(temp)
        return temp[ip], kouho[ip]
    end

    function fmin()
        kouho = lx[1:ni]
        rdiag = [r[i, i] for i in kouho]
        temp = r[q1, kouho] .^ 2 ./ (rdiag .* r[q1, q1])
        ip = argmin(temp)
        return temp[ip], lx[ip]
    end

    function sweepsreg!(r, ip, q1)
        ap = r[ip, ip]
        abs(ap) > EPSINV || error("正規方程式の係数行列が特異行列です")
        for i in 1:q1
            if i != ip
                temp = r[ip, i] / ap
                for j in 1:q1
                    if j != ip
                        r[j, i] -= r[j, ip] * temp
                    end
                end
            end
        end
        r[:, ip] /= ap
        r[ip, :] /= -ap
        r[ip, ip] = 1 / ap
    end

    function printsreg()
        variable = results[:variables]
        @printf("\n%-11s  %10s  %10s  %9s  %9s  %10s  %9s\n",
                 "", "Estimate", "Std. Error", "t value", "Pr(>|t|)", "Std. coef.", "Tolerance")
        for i in 1:length(variable)
            @printf("%-11s  %10.6f  %10.6f  %9.6f  %9.6f  %10.6f  %9.6f\n",
                    variable[i], results[:coefficient][i], results[:standarderror][i],
                    results[:tvalue][i], results[:tPvalue][i],
                    results[:standardizedcoefficient][i],
                    results[:tolerance][i])
        end
        @printf("%-11s  %10.6f  %10.6f  %9.6f  %9.6f\n\n",
                "(Intercept)", results[:coefficient][end], results[:standarderror][end],
                results[:tvalue][end], results[:tPvalue][end])
    end

    function additional()
        @printf("Residual standard error: %.5f on %d degrees of freedom\n",
                sqrt(results[:ve]), results[:df][2])
        @printf("Multiple correlation coefficient: %.5f\n",
                results[:R])
        @printf("Multiple R-squared: %.5f, Adjusted R-squared: %.5f\n",
                results[:R2], results[:adjR2])
        @printf("log likelifood: %.5f,  AIC: %.5f\n",
                results[:loglik], results[:AIC])
    end

    function printanova()
        variable = results[:variables]
        @printf("\nANOVA Table\n%11s  %10s %4s  %9s  %10s  %8s\n",
                "", "SS", "df", "MS", "F value", "p value")
        @printf("%-11s  %10.6f %4d  %9.6f  %10.6f  %8.6f\n",
                "Regression", results[:SS][1], results[:df][1], results[:MS][1], results[:Fvalue], results[:FPvalue])
        @printf("%-11s  %10.6f %4d  %9.6f\n",
                "Residual", results[:SS][2], results[:df][2], results[:MS][2])
        @printf("%-11s  %10.6f %4d\n",
                "Total", results[:SS][3], results[:df][3])
    end

    function procpredict(predict)
        lxi = results[:lxi]
        ve = results[:ve]
        results[:coefficient]
        constant = results[:coefficient][end]
        b = results[:coefficient][1:end-1]
        y = data[:, size(data, 2)]
        data = data[:, lxi]
        est = data * b .+ constant
        res = y .- est
        r = r[lxi, lxi]
        means = results[:mean][lxi]
        stds = results[:sd][lxi]
        nc, q = size(data)
        factor = sqrt(nc / (nc - 1))
        for i in 1:q
            data[:, i] = ((data[:, i] .- means[i]) ./ stds[i]) .* factor
        end
        stres = [transpose(data[i, :]) * r * data[i, :] for i = 1:nc]
        stres = res ./ sqrt.(ve .* (1 .- (stres .+ 1) ./ nc))
        if predict
            @printf("\npredict\n%5s %12s  %12s  %12s  %12s\n",
                    "i", "observed", "predicted", "residuals", "std. resid.")
            for i = 1:nc
                @printf("%5d:%12.6f  %12.6f  %12.6f  %12.6f\n",
                        i, y[i], est[i], res[i], stres[i])
            end
        end
        results = merge(results, Dict(:observed => y, :predicted => est,
                    :residuals => res, :standardizedresiduals => stres))
    end

    EPSINV = 1e-6
    MULTICO = 10
    Pout >= Pin || (Pout = Pin)
    step = 0
    oldrhoy = 1
    ve = 0
    ip = 0
    lxi = []
    nc, q1 = size(data)
    q = q1 - 1
    (stepwise == true || nc > q1) || error("独立変数の個数よりデータ数が多くなければなりません")
    q > 1 || (stepwise = false)
    length(vname) != 0 || (vname = vcat(["x" * string(i) for i = 1:q], "y"))
    means = vec(mean(data, dims=1));
    sd = vec(std(data, dims=1));
    var = sd .^ 2;
    ss = (cov(data) .* (nc - 1))[:, q1];
    r = cor(data);
    results = Dict(:samplesize => nc, :dependentvariable => vname[q1],
                   :independentvariable => vname[1:q],
                   :mean => means, :var => var, :sd => sd, :r => r)
    stde = similar(means);
    tv = similar(means);
    pp = similar(means);
    ni = 0
    if stepwise == false
        for ip in 1:q
            sweepsreg!(r, ip, q1)
        end
        lx = 1:q
        ni = q
        stepout(0)
    else
        if verbose
            println("\n変数編入基準    Pin:  $Pin")
            println("変数除去基準    Pout: $Pout\n")
        end
        lx = zeros(Int, q)
        ni = 0
        while ni < min(q, nc - 2)
            ansmax, ip = fmax()
            P = (nc - ni - 2) * ansmax # [1]
            P = pf(P, 1, nc - ni - 2, false)
            verbose && print("編入候補変数: $(vname[ip])   P : $(formatpval(P))")
            if P > Pin
                verbose && println("  編入されませんでした")
                break
            end
            verbose && println("  編入されました")
            ni = ni + 1
            lx[ni] = ip
            sweepsreg!(r, ip, q1)
            stepout(1)
            printsreg()
            while true
                ansmin, ip = fmin()
                P = (nc - ni - 1) * ansmin # [1]
                P = pf(P, 1, nc - ni - 1, false)
                verbose && print("除去候補変数: $(vname[ip])   P :$(formatpval(P))")
                if P <= Pout
                    verbose && println("  除去されませんでした")
                    break
                else
                    verbose && println("  除去されました")
                    lx = lx[ - which(lx == ip)]
                    ni = ni - 1
                    sweepsreg!(r, ip, q1)
                    stepout(2)
                    printsreg()
                end
            end
        end
    end

    ni == 0 && error("条件（ Pin < $Pin）を満たす独立変数がありません")
    println("\n***** 最終結果 *****")
    printsreg()
    additional()
    printanova()
    procpredict(predict)
end

using RDatasets
iris = dataset("datasets", "iris");
a = iris[:, 1:4];
result = sreg(Array(a), stepwise=true,
    vname=["SepalLength", "SepalWidth", "PetalLength", "PetalWidth"],
    predict=true, verbose=true)

#==========
    変数編入基準    Pin:  0.05
    変数除去基準    Pout: 0.05

    編入候補変数: PetalLength   P : < 0.0001  編入されました
    1 lxi = [3]

     ***** ステップ 1 *****   　　編入変数: PetalLength

                   Estimate  Std. Error    t value   Pr(>|t|)  Std. coef.  Tolerance
    PetalLength    0.415755    0.009582  43.387237   0.000000    0.962865   1.000000
    (Intercept)   -0.363076    0.039762   9.131221   0.000000

    除去候補変数: PetalLength   P :< 0.0001  除去されませんでした
    編入候補変数: SepalWidth   P : 0.02014  編入されました
    1 lxi = [3, 2]

     ***** ステップ 2 *****   　　編入変数: SepalWidth

                   Estimate  Std. Error    t value   Pr(>|t|)  Std. coef.  Tolerance
    PetalLength    0.426270    0.010447  40.803902   0.000000    0.987217   0.816439
    SepalWidth     0.099396    0.042310   2.349216   0.020144    0.056837   0.816439
    (Intercept)   -0.706478    0.151334   4.668332   0.000007

    除去候補変数: SepalWidth   P :0.02014  除去されませんでした
    編入候補変数: SepalLength   P : < 0.0001  編入されました
    1 lxi = [3, 2, 1]

     ***** ステップ 3 *****   　　編入変数: SepalLength

                   Estimate  Std. Error    t value   Pr(>|t|)  Std. coef.  Tolerance
    PetalLength    0.524083    0.024491  21.398708   0.000000    1.213746   0.132314
    SepalWidth     0.222829    0.048938   4.553279   0.000011    0.127419   0.543585
    SepalLength   -0.207266    0.047506   4.362929   0.000024   -0.225166   0.159822
    (Intercept)   -0.240307    0.178370   1.347243   0.179989

    除去候補変数: SepalLength   P :< 0.0001  除去されませんでした

    ***** 最終結果 *****

                   Estimate  Std. Error    t value   Pr(>|t|)  Std. coef.  Tolerance
    PetalLength    0.524083    0.024491  21.398708   0.000000    1.213746   0.132314
    SepalWidth     0.222829    0.048938   4.553279   0.000011    0.127419   0.543585
    SepalLength   -0.207266    0.047506   4.362929   0.000024   -0.225166   0.159822
    (Intercept)   -0.240307    0.178370   1.347243   0.179989

    Residual standard error: 0.19197 on 146 degrees of freedom
    Multiple correlation coefficient: 0.96843
    Multiple R-squared: 0.93785, Adjusted R-squared: 0.93657
    log likelifood: 36.75109,  AIC: -63.50218

    ANOVA Table
                         SS   df         MS     F value   p value
    Regression    81.189636    3  27.063212  734.388537  0.000000
    Residual       5.380298  146   0.036851
    Total         86.569933  149

    predict
        i     observed     predicted     residuals   std. resid.
        1:    0.200000      0.216252     -0.016252     -0.085561
        2:    0.200000      0.146291      0.053709      0.283384
        3:    0.200000      0.179901      0.020099      0.105791
        4:    0.200000      0.283162     -0.083162     -0.438018
        5:    0.200000      0.259261     -0.059261     -0.312245
        6:    0.400000      0.400428     -0.000428     -0.002269
        7:    0.300000      0.297602      0.002398      0.012656
        8:    0.200000      0.267104     -0.067104     -0.352773
        9:    0.200000      0.227641     -0.027641     -0.146141
       10:    0.100000      0.220982     -0.120982     -0.636584
                :
==========#

頑張ってます。クリックお願いします。

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

Julia に翻訳--145　重回帰分析

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

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

重回帰分析
http://aoki2.si.gunma-u.ac.jp/R/mreg.html

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

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

出力にこだわらなければ簡単なプログラムなんだけど。

==========#

using Statistics, LinearAlgebra, Rmath, Printf

function mreg(dat)
    n, nc = size(dat)
    nv = nc - 1
    r = cor(dat)
    mall = vec(mean(dat, dims=1))
    m = mall[1:nc-1]
    betas = r[1:nc-1, 1:nc-1] \ r[1:nc-1, nc]
    variance = cov(dat) .* (n - 1)
    prop = diag(variance)
    prop = (prop / prop[nc])[1:nc-1]
    b = betas ./ sqrt.(prop)
    Sr = transpose(variance[:, nc][1:nc-1]) * b
    St = variance[nc, nc]
    Se = St .- Sr
    SS = [Sr, Se, St]
    dfr = nv
    dfe = n - nv - 1
    dft = n - 1
    df = [dfr, dfe, dft]
    Ms = SS ./ df
    f = Ms[1] / Ms[2]
    fvalue = [f, NaN, NaN]
    p = [pf(f, dfr, dfe, false),NaN, NaN]
    b0 = mall[nc] - sum(b .* m)
    b = vcat(b, b0)
    inverse = inv((n - 1) .* cov(dat)[1:nc-1, 1:nc-1])
    SEb = vcat([sqrt(inverse[i, i] * Ms[2]) for i =1:nv],
               sqrt((1 / n + transpose(m) * inverse * m) * Ms[2]))
    tval = b ./ SEb
    pval = pt.(abs.(tval), n - nv - 1, false) * 2
    tolerance = 1 ./ diag(inv(cor(dat)[1:nc-1, 1:nc-1]))
    R2 = 1 - Se / St
    R = sqrt(R2)
    R2s = 1 - Ms[2] / Ms[3]
    loglik = - 0.5 * n * (log(2 * pi) + 1 - log(n) + log(Se))
    AIC = 2 * nc + 2 - 2 * loglik
    printresult(b, SEb, tval, pval, betas, tolerance)
    printanovatable(SS, df, Ms, f, p)
    printstatistics(R2, R, R2s, loglik, AIC)
    Dict(:result => hcat(b, SEb, tval, pval, vcat(betas, NaN), vcat(tolerance, NaN)),
            :anova => hcat(SS, df, Ms, fvalue, p),
            :Rs => [R2, R, R2s, loglik, AIC])
end

diag(a) = [a[i, i] for i in 1:size(a, 1)]

function printresult(b, SEb, tval, pval, betas, tolerance)
    @printf("%6s%12s%12s%12s%12s%18s%12s\n",
            "", "偏回帰係数", "標準誤差", "t 値", "P 値", "標準化偏回帰係数", "トレランス")
    for i = 1:length(betas)
        @printf("%6s%12.5g%12.5g%12.5g%12.5g%18.5g%12.5g\n",
                "Var" * string(i), b[i], SEb[i], tval[i], pval[i], betas[i], tolerance[i])
    end
    @printf("%6s%12.5g%12.5g%12.5g%12.5g\n",
            "定数項", b[end], SEb[end], tval[end], pval[end])
end

function printanovatable(SS, df, Ms, f, p)
    println("回帰の分散分析表")
    @printf("%4s  %10s  %6s  %12s  %12s  %12s\n",
               "", "平方和", "自由度", "平均平方", "F 値", "P 値")
    @printf("回帰  %10.6g  %6d  %12.6g  %12.6g  %12.6g\n",
               SS[1], df[1], Ms[1], f[1], p[1])
    @printf("残差  %10.6g  %6d  %12.6g\n", SS[2], df[2], Ms[2])
    @printf("全体  %10.6g  %6d  %12.6g\n", SS[3], df[3], Ms[3])
end

function printstatistics(R2, R, R2s, loglik, AIC)
    println("重相関係数 = $R")
    println("重相関係数の二乗 = $R2")
    println("自由度調整済重相関係数の二乗 = $R2s")
    println("対数尤度 = $loglik")
    println("AIC = $AIC")
end

dat = [
    1.2 1.9 0.9
    1.6 2.7 1.3
    3.5 3.7 2.0
    4.0 3.1 1.8
    5.6 3.5 2.2
    5.7 7.5 3.5
    6.7 1.2 1.9
    7.5 3.7 2.7
    8.5 0.6 2.1
    9.7 5.1 3.6]

mreg(dat)
# 偏回帰係数    標準誤差        t 値        P 値  標準化偏回帰係数  トレランス
# Var1     0.20462   0.0075643      27.051  2.4192e-08           0.67067     0.98358
# Var2     0.28663    0.010802      26.536  2.7638e-08           0.65793     0.98358
# 定数項     0.14918    0.054506      2.7368     0.02905
# 回帰の分散分析表
#   平方和  自由度      平均平方          F 値          P 値
# 回帰     6.67164       2       3.33582       823.477   4.93187e-09
# 残差   0.0283563       7     0.0040509
# 全体         6.7       9      0.744444
# 重相関係数 = 0.9978816139144453
# 重相関係数の二乗 = 0.9957677153884982
# 自由度調整済重相関係数の二乗 = 0.9945584912137834
# 対数尤度 = 15.13806917315012
# AIC = -22.27613834630024

頑張ってます。クリックお願いします。

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

Julia に翻訳--144　リッカート尺度

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

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

リッカート尺度
http://aoki2.si.gunma-u.ac.jp/R/likert.html

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

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

==========#

using Rmath

function likert(dat)
    resp = dat / sum(dat)
    cum = cumsum(resp)
    result =(dnorm.(qnorm.(vcat(0, cum[1:end-1]))) .- dnorm.(qnorm.(cum))) ./ resp
    println(result)
    Dict(:result => result)
end

dat = [7, 18, 34, 26, 15]
likert(dat)
# [-1.9181130652891667, -1.0194925450770291, -0.20873614063905502, 0.5984157124911674, 1.5543918350245474]

頑張ってます。クリックお願いします。

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

Julia に翻訳--143　クロンバックのα信頼性係数

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

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

クロンバックのα信頼性係数
http://aoki2.si.gunma-u.ac.jp/R/alpha.html

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

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

出力にはこだわらない。
必要なら，戻り値を使って。

==========#

using Statistics

function alpha(x; detail = false)
    nr, k = size(x)
    (k > 1 && nr > 1) || error("parameter error")
    α = alpha0(x)
    println("α = $α")
    if detail == true && k >= 3
        R2 = 1 .- 1 ./ diag(inv(cor(x)))
        α2 = similar(R2)
        cor2 = similar(R2)
        for i in 1:k
            x2 = x[:, [i != j for j in 1:k]]
            α2[i] = alpha0(x2)
            cor2[i] = cor(x[:, i], vec(sum(x2, dims=2)))
        end
        println("それぞれの変数を除いたときの α\n$α2")
        println("それぞれの変数とその変数を除いたときの合計値との相関係数\n$cor2")
        println("それぞれの変数をその変数以外の変数で予測したときの決定係数\n$R2")
        Dict(:alpha => α, :alpha2 => α2, :R2 => R2, :cor2 => cor2)
    end
end

function alpha0(x)
    k = size(x, 2)
    VarCovMat = cov(x)
    Sy2 = sum(VarCovMat)
    Sj2 = sum(diag(VarCovMat))
    k / (k - 1) * (1 - Sj2 / Sy2)
end

diag(a) = [a[i, i] for i in 1:size(a, 1)]

x = [
    49 44 37 54
    36 36 36 42
    42 51 45 35
    47 67 54 40
    54 59 68 54
    51 55 59 67
    45 36 48 46
    72 49 50 58];

alpha(x)
alpha(x, detail = true)
# α = 0.7403365794138637
# それぞれの変数を除いたときの α
# [0.6657399977072106, 0.7307558528428093, 0.592778807780536, 0.72184]
# それぞれの変数とその変数を除いたときの合計値との相関係数
# [0.5599374682949105, 0.4444039960775981, 0.6801937600489253, 0.4589177292647371]
# それぞれの変数をその変数以外の変数で予測したときの決定係数
# [0.4298971517659005, 0.551857280000196, 0.6071655584759322, 0.4821633678144336]

頑張ってます。クリックお願いします。

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

Julia に翻訳--142　スミルノフ・グラブス検定

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

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

スミルノフ・グラブス検定
http://aoki2.si.gunma-u.ac.jp/R/SG.html

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

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

==========#

using Statistics, Rmath

function sg(x)
    n = length(x)
    t = abs.(extrema(x) .- mean(x)) ./ std(x)
    df = n - 2
    p = @. n * pt(sqrt(df / ((n - 1) ^ 2 / t ^ 2 / n - 1)), df, false)
    p = [min(p0, 1) for p0 in p]
    println("min: $(minimum(x))")
    println("t = $(t[1]),  df = $df,  pvalue = $(p[1])")
    println("max: $(maximum(x))")
    println("t = $(t[2]),  df = $df,  pvalue = $(p[2])")
    Dict(:t => t, :df => df, :p => p)
end

x = [133, 134, 134, 134, 135, 135, 139, 140, 140, 140,
    141, 142, 142, 144, 144, 147, 147, 149, 150, 164]

sg(x)
# min: 133
# t = 1.1724347855313042,  df = 18,  pvalue = 1.0
# max: 164
# t = 3.0052064042928888,  df = 18,  pvalue = 0.004864034251605727

 

頑張ってます。クリックお願いします。

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

Julia に翻訳--141　Breslow-Day 検定

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

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

Breslow-Day 検定
http://aoki2.si.gunma-u.ac.jp/R/BD-test.html

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

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

@. なんか書かなくてもよいようにしてほしいなあ

==========#

using Rmath

function bdtest(m)
    d1, d2, k = size(m)
    Nk = vec(sum(m, dims=[1,2])) # vex() がないと，1x1x8 の3次元行列になる
    psiMH = sum(m[1, 1, :] .* m[2, 2, :] ./ Nk) / sum(m[2, 1, :] .* m[1, 2, :] ./ Nk)
    nk1 = m[1, 1, :] .+ m[1, 2, :]
    nk2 = m[2, 1, :] .+ m[2, 2, :]
    xkt = m[1, 1, :] .+ m[2, 1, :]
    a1 = psiMH - 1
    @. b1 = -psiMH * (nk1 + xkt) - nk2 + xkt
    @. c1 = psiMH * nk1 * xkt
    @. e = (-b1 - sqrt(b1 ^ 2 - 4 * a1 * c1)) / (2 * a1)
    @. v = 1 / (1 / e + 1 / (nk1 - e) + 1 / (xkt - e) + 1 / (nk2 - xkt + e))
    chisqBD = sum((m[1, 1, :] .- e) .^ 2 ./ v) # この行は @. を使うとエラーになる
    df = k - 1
    p = pchisq(chisqBD, df, false)
    println("chisq. = $chisqBD,  df = $df,  p value = $p")
    Dict(:chisq => chisqBD, :df => df, :pvalue => p)
end

m = [1, 8, 0, 10,
     2, 47, 0, 60,
     3, 29, 1, 13,
     3, 17, 0, 7,
     13, 45, 1, 45,
     13, 26, 0, 18,
     8, 15, 0, 10,
     11, 4, 0, 4 ];
m = reshape(m, 2, 2, :)

bdtest(m)
# chisq. = 8.943202421541203,  df = 7,  p value = 0.25675965693687036

 

頑張ってます。クリックお願いします。

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

Julia の小ネタ--014　三次元配列の作り方

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

以下のような 2×3×4 行列を作る方法は 2 通りある。

2×3×4 Array{Int64,3}:
[:, :, 1] =
 1  3  5
 2  4  6

[:, :, 2] =
 7   9  11
 8  10  12

[:, :, 3] =
 13  15  17
 14  16  18

[:, :, 4] =
 19  21  23
 20  22  24

第1の方法

Julia は R や FORTRAN と同じように，列優先（次元表示の左にある添え字が最も速く変化する）なので，列優先で 1 次元配列（ベクトル）を用意して，reshape() で次元を定義し直す。

reshape([1,2,3,4,5,6,7,8,9,10,
         11,12,13,14,15,16,17,18,19,20,
         21,22,23,24], (2, 3, 4))

この場合は，以下のように書くこともできる。

reshape(1:24, 2, 3, :)

第2の方法

行優先で第1，第2次元の配列を定義し，cat で連結するときに dims=3 を指定する。

cat([1 3 5; 2 4 6], [7 9 11; 8 10 12],
    [13 15 17; 14 16 18], [19 21 23; 20 22 24], dims=3)

 

頑張ってます。クリックお願いします。

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

Julia に翻訳--140　Chow 検定

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

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

Chow 検定
http://aoki2.si.gunma-u.ac.jp/R/Chow.html

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

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

==========#

using Statistics, Rmath

chow = function(dat1, dat2)
    ess12 = ess(dat1) + ess(dat2)
    essc = ess(vcat(dat1, dat2))
    nr, df1 = size(dat1)
    df2 = nr + size(dat2, 1) - 2 * df1
    f = (essc - ess12) * df2 / (df1 * ess12)
    p = pf(f, df1, df2, false)
    println("F = $f,  df1 = $df1,  df2 = $df2,  p value = $p")
    Dict(:F => f, :df1 => df1, :df2 => df2, :pvalue => p)
end

function ess(dat) # 回帰分析を行い，誤差変動を返す関数
    r = cor(dat)
    betas = r[1:end-1, 1:end-1] \ r[1:end-1, end]
    SS = cov(dat, corrected=false) .* size(dat, 1)
    prop = [SS[i, i] for i in 1:size(SS, 1)]
    Sr = sum(SS[1:end-1, end] .* (betas .* sqrt.(prop[end] ./ prop[1:end-1])))
    St = SS[end, end]
    St - Sr
end

dat1 = [
    1.2 1.9 0.9
    1.6 2.7 1.3
    3.5 3.7 2.0
    4.0 3.1 1.8
    5.6 3.5 2.2
    5.7 7.5 3.5
    6.7 1.2 1.9
    7.5 3.7 2.7
    8.5 0.6 2.1
    9.7 5.1 3.6
    ];

dat2 = [
    1.4 1.3 0.5
    1.5 2.3 1.3
    3.1 3.2 2.5
    4.4 3.6 1.1
    5.1 3.1 2.8
    5.2 7.3 3.3
    6.5 1.5 1.3
    7.8 3.2 2.2
    8.1 0.1 2.8
    9.5 5.6 3.9
    ];

chow(dat1, dat2)
# F = 0.07154752544573709,  df1 = 3,  df2 = 14,  p value = 0.9742318542171906

頑張ってます。クリックお願いします。

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

Julia に翻訳--139　分散・共分散行列の同等性の検定

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

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

分散・共分散行列の同等性の検定
http://aoki2.si.gunma-u.ac.jp/R/eq-cov.html

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

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

関数を入れ子にして一行で書くのもよいが，for ループで書くとわかりやすいというメリットもある。

==========#

using Statistics, LinearAlgebra, Rmath

function boxm(x, gvar)
    nv = size(x, 2)
    nv >= 2 || error("分散共分散行列を計算する変数は2個以上必要です")
    index, ni = table(gvar)
    n = length(gvar)
    g = length(ni)
    sumlogdetSi = 0
    S = zeros(nv, nv)
    for (i, gvari) in enumerate(index)
        y = x[gvar .== gvari, :]
        covy = cov(y)
        sumlogdetSi += (ni[i] - 1) * log(det(covy))
        S .+= (ni[i] - 1) .* covy
    end
    S ./= n - g
    M = (n - g) * log(det(S)) - sumlogdetSi
    df1 = (g - 1) * nv * (nv + 1) / 2
    rho = 1 - (2 * nv ^ 2 + 3 * nv - 1) / (6 * (nv + 1) * (g - 1)) * (sum(1 ./ (ni .- 1)) - 1 / (n - g))
    tau = (nv - 1) * (nv + 2) / (6 * (g - 1)) * (sum(1 ./ (ni .- 1) .^ 2) - 1 / (n - g) ^ 2)
    df2 = (df1 + 2) / abs(tau - (1 - rho) ^ 2)
    gamma = (rho - df1 / df2) / df1
    F = M * gamma
    p = pf(F, df1, df2, false)
    println("M = $M\nF = $F,  df1 = $df1,  df2 = $df2\np value = $p")
    Dict(:M => M, :F => F, :df1 => df1, :df2 => df2, :pvalue => p)
end

function table(x) # indices が少ないとき
    indices = sort(unique(x))
    counts = zeros(Int, length(indices))
    for i in indexin(x, indices)
        counts[i] += 1
    end
    return indices, counts
end

function rep(x, n::Array{Int64,1})
    length(x) == length(n) || error("length(x) wasn't length(n)")
    vcat([repeat([x[i]], n[i]) for i = 1:length(x)]...)
end

x = [
   2.9 161.7 120.8
   2.3 114.8  85.2
   2.0 128.4  92.0
   3.2 149.2  97.3
   2.7 126.0  81.1
   4.4 133.8 107.6
   4.1 161.3 114.0
   2.1 111.5  77.3
   4.8 198.7 172.9
   3.6 199.3 157.9
   2.0 188.4 152.7
   4.9 183.6 164.2
   3.9 173.5 172.2
   4.4 184.9 163.2 ];
gvar = rep(1:2, [8, 6]);

boxm(x, gvar)
# M = 9.730421574785467
# F = 1.1535477771673381,  df1 = 6.0,  df2 = 795.1946969489453
# p value = 0.32937843690375435

頑張ってます。クリックお願いします。

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

Julia に翻訳--138　Brown-Forsythe 検定，三群以上の場合の等分散性の検定

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

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

Brown-Forsythe 検定（三群以上の場合の等分散性の検定）
http://aoki2.si.gunma-u.ac.jp/R/Brown-Forsythe-test.html

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

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

R の split() と sapply() の間単版を Julia で書いてみた。

==========#

using Statistics, Rmath

function brownforsythetest(x, group)
    d = split2(x, group)
    df1 = length(d) - 1
    ni = sapply(d, length)
    mi = sapply(d, mean)
    ui = sapply(d, var)
    ci = (1 .- ni ./ sum(ni)) .* ui
    FBF = sum(ni .* (mi .- mean(x)) .^ 2) / sum(ci)
    C = ci ./ sum(ci)
    df2 = 1 / sum(C .^ 2 / (ni .- 1))
    p = pf(FBF, df1, df2, false)
    println("F = $FBF,  df1 = $df1,  df2 = $df2,  p value = $p")
    Dict(:F => FBF, :df1 => df1, :df2 => df2, :pvalue => p)
end

function split2(x, g)
    index = sort(unique(g))
    y = Array{Array{Float64,1},1}(undef, length(index))
    [y[i] = x[g .== i] for i in index]
end

sapply(x, FUNC) = map(FUNC, x)

function rep(x; each::Int64)
    vcat([repeat([x[i]], each) for i = 1:length(x)]...)
end

x = [3, 3, 4, 2, 5, 2, 3, 4, 8, 8, 5, 6];
g = rep(1:3, each = 4);
brownforsythetest(x, g)
# F = 10.854545454545452,  df1 = 2,  df2 = 7.606873428331936,  p value = 0.00590981731767559

頑張ってます。クリックお願いします。

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

Julia に翻訳--137　Levene 検定，三群以上の場合の等分散性の検定

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

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

Levene 検定（三群以上の場合の等分散性の検定）
http://aoki2.si.gunma-u.ac.jp/R/levene-test.html

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

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

==========#

using Statistics, Rmath

function levenetest(x, group; method = "mean") # or method = "median"
    n = length(x)
    ni = tapply(x, group, length)
    k = length(ni)
    FUNC = method == "mean" ? mean : median
    x = abs.(x .- tapply(x, group, FUNC)[group])
    sw = sum((ni .- 1) .* tapply(x, group, var))
    dfw = n - k
    dfb = k - 1
    f = ((var(x) * (n - 1) - sw) / dfb) / (sw / dfw)
    P = pf(f, dfb, dfw, false)
    println("等分散性の検定（$method で調整）")
    println("F = $f,  dfb = $dfb,  dfw = $dfw,  p value = $P")
    Dict(:F => f, :dfb => dfb, :dfw => dfw, :pvalue => P, :method => method)
end

function tapply(x, g, func)
    indices = sort(unique(g))
    [func(x[g .== i]) for i in indices]
end

function rep(x; each::Int64)
    vcat([repeat([x[i]], each) for i = 1:length(x)]...)
end

x = [3, 3, 4, 2, 5, 2, 3, 4, 8, 8, 5, 6]
g = rep(1:3, each = 4)

levenetest(x, g)
# 等分散性の検定（mean で調整）
# F = 2.0999999999999996,  dfb = 2,  dfw = 9,  p value = 0.17844673330519398

levenetest(x, g, method = "median")
# 等分散性の検定（median で調整）
# F = 1.9090909090909094,  dfb = 2,  dfw = 9,  p value = 0.203644351599495

頑張ってます。クリックお願いします。

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

Julia に翻訳--136　一元配置分散分析，exact test

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

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

一元配置分散分析(exact test)
http://aoki2.si.gunma-u.ac.jp/R/exact-oneway-test.html

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

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

==========#

using Statistics, Rmath, Printf

function exactonewaytest(x, y=[]; permutation = true)
    function found()
        hh = performtest(um)
        if hh <= pvalue + EPSILON
            nprod = sum(perm_table[rt .+ 1]) - sum(perm_table[um .+ 1])
            nntrue += exp(nprod - nntrue2 * log_expmax)
            while nntrue >= EXPMAX
                nntrue /= EXPMAX
                nntrue2 += 1
            end
        end
        ntables += 1
    end

    function search(x, y)
        if y == 1
            found()
        elseif x == 1
            t = um[1, 1] - um[y, 1]
            if t >= 0
                um[1, 1] = t
                search(nc, y - 1)
                um[1, 1] += um[y, 1]
            end
        else
            search(x - 1, y)
            while um[y, 1] != 0 && um[1, x] != 0
                um[y, x] += 1
                um[y, 1] -= 1
                um[1, x] -= 1
                search(x - 1, y)
            end
            um[y, 1] += um[y, x]
            um[1, x] += um[y, x]
            um[y, x] = 0
        end
    end

    function permutationtest()
        denom2 = 0
        denom = perm_table[n + 1] - sum(perm_table[ct .+ 1])
        while denom > log_expmax
            denom -= log_expmax
            denom2 += 1
        end
        denom = exp(denom)
        um[:, 1] = rt
        um[1, :] = ct
        search(nc, nr)
        pvalue = nntrue / denom * EXPMAX ^ (nntrue2 - denom2)
        @printf("並べ替え検定による P 値 = % .10g\n", pvalue)
        @printf("査察した分割表の個数は % s 個\n", ntables)
        return pvalue
    end

    function performtest(u)
        x = Array{Array{Float64,1},1}(undef, nr)
        for i in 1:nr
            x[i] = score[u[i, :] .> 0]
        end
        onewaytest(x)
    end

    function simpleonewaytest()
        @printf("観察値による一元配置分散分析の P 値 = % g\n", pvalue)
    end

    if length(y) == 0
        ni = map(length, x)
        y = rep(1:length(ni), ni)
        indexy, score, t = table(y, vcat(x...))
    else
        indexy, score, t = table(y, x)
    end
    EPSILON = 1e-10
    EXPMAX = 1e100
    log_expmax = log(EXPMAX)
    nr, nc = size(t)
    um = zeros(Int, nr, nc)
    rt = sum(t, dims=2)
    ct = transpose(sum(t, dims=1))
    n = sum(t)
    q = cumsum(vcat(0, ct[1:nc-1])) .+ (ct .+ 1) .* 0.5
    half = (n + 1) * 0.5
    pvalue = performtest(t)
    simpleonewaytest()
    if permutation
        perm_table = cumsum(vcat(0, log.(1:n + 1)))
        ntables = nntrue = nntrue2 = 0
        permutationtest()
    end
end

function onewaytest(x)
    n = map(length, x)
    m = map(mean, x)
    v = map(var, x)
    w = n ./ v
    sumw = sum(w)
    tmp = sum((1 .- (w / sumw)) .^ 2 ./ (n .- 1)) / (nr ^ 2 - 1)
    m0 = sum(w .* m) / sumw
    F = sum(w .* (m .- m0) .^ 2) / ((nr - 1) * (1 + 2 * (nr - 2) * tmp))
    df1 = nr - 1
    df2 = 1 / (3 * tmp)
    pf(F, df1, df2, false)
end

function rep(x, n::Array{Int64,1})
    length(x) == length(n) || error("length(x) wasn't length(n)")
    vcat([repeat([x[i]], n[i]) for i = 1:length(x)]...)
end

function table(x, y) # 二次元
    indicesx = sort(unique(x))
    indicesy = sort(unique(y))
    counts = zeros(Int, length(indicesx), length(indicesy))
    for (i, j) in zip(indexin(x, indicesx), indexin(y, indicesy))
        counts[i, j] += 1
    end
    return indicesx, indicesy, counts
end

x = [[36.7, 52.4, 65.8], [45.7, 61.9, 65.3], [52.6, 76.6, 81.3]];
exactonewaytest(x)
# 観察値による一元配置分散分析の P 値 =  0.440037
# 並べ替え検定による P 値 =  0.4107142857
# 査察した分割表の個数は 1680 個

 

頑張ってます。クリックお願いします。

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

Julia に翻訳--135　クラスカル・ウォリス検定，exact test

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

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

クラスカル・ウォリス検定(exact test)
http://aoki2.si.gunma-u.ac.jp/R/exact-kw.html

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

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

global は不要

==========#

using Rmath, Printf

function exactkw(x, y = []; exact = true)
    function found()
        hh = sum((um * q) .^ 2 ./ rt)
        if hh >= stat_val || abs(hh - stat_val) <= EPSILON
            nprod = sum(perm_table[rt .+ 1]) - sum(perm_table[um .+ 1])
            nntrue += exp(nprod - nntrue2 * log_expmax)
            while nntrue >= EXPMAX
                nntrue /= EXPMAX
                nntrue2 += 1
            end
        end
        ntables += 1
    end

    function search(x, y)
        if y == 1
            found()
        elseif x == 1
            t = um[1, 1] - um[y, 1]
            if t >= 0
                um[1, 1] = t
                search(nc, y - 1)
                um[1, 1] += um[y, 1]
            end
        else
            search(x - 1, y)
            while um[y, 1] != 0 && um[1, x] != 0
                um[y, x] += 1
                um[y, 1] -= 1
                um[1, x] -= 1
                search(x - 1, y)
            end
            um[y, 1] += um[y, x]
            um[1, x] += um[y, x]
            um[y, x] = 0
        end
    end

    function exacttest()
        denom2 = 0
        denom = perm_table[n + 1] - sum(perm_table[ct .+ 1])
        while denom > log_expmax
            denom -= log_expmax
            denom2 += 1
        end
        denom = exp(denom)
        um[:, 1] = rt
        um[1, :] = ct
        search(nc, nr)
        @printf("正確な P 値 = % .10g\n", nntrue / denom * EXPMAX ^ (nntrue2 - denom2))
        @printf("査察した分割表の個数は % s 個\n", ntables)
    end

    function kwtest(u)
        return sum((u * q) .^ 2 ./ rt)
    end

    function asymptotic()
        chisq = (stat_val * 12 / (n * (n + 1)) - 3 * (n + 1)) / (1 - sum(ct .^ 3 .- ct) / (n ^ 3 - n))
        @printf("カイ二乗値 = % g, 自由度 = % i, P 値 = % g\n", chisq, nr - 1, pchisq(chisq, nr - 1, false))
    end

    if ndims(x) == 2
        t = x
    elseif length(y) == 0
        ni = map(length, x)
        y = rep(1:length(ni), ni)
        indexy, indexx, t = table(y, vcat(x...))
    else
        indexy, indexx, t = table(y, x)
    end
    EPSILON = 1e-10
    EXPMAX = 1e100
    log_expmax = log(EXPMAX)
    nr, nc = size(t)
    rt = sum(t, dims=2)
    ct = transpose(sum(t, dims=1))
    um = zeros(Int, nr, nc)
    n = sum(t)
    q = cumsum(vcat(0, ct[1:nc-1])) .+ (ct .+ 1) .* 0.5
    half = (n + 1) * 0.5
    stat_val = kwtest(t)
    asymptotic()
    if exact
        perm_table = cumsum(vcat(0, log.(1:(n + 1))))
        ntables = nntrue = nntrue2 = 0
        exacttest()
    end
end

function rep(x, n::Array{Int64,1})
    length(x) == length(n) || error("length(x) wasn't length(n)")
    vcat([repeat([x[i]], n[i]) for i = 1:length(x)]...)
end

function table(x, y) # 二次元
    indicesx = sort(unique(x))
    indicesy = sort(unique(y))
    counts = zeros(Int, length(indicesx), length(indicesy))
    for (i, j) in zip(indexin(x, indicesx), indexin(y, indicesy))
        counts[i, j] += 1
    end
    return indicesx, indicesy, counts
end

分割表を与える場合 Array{Int64,2}

x = [5 3 2 1
     4 3 5 2
     2 3 1 2]
exactkw(x)
# カイ二乗値 =  1.32485, 自由度 =  2, P 値 =  0.5156
# 正確な P 値 =  0.5268191237
# 査察した分割表の個数は 24871 個

データベクトルと factor ベクトルを与える場合

data = [
    3.42, 3.84, 3.96, 3.76,
    3.17, 3.63, 3.47, 3.44, 3.39,
    3.64, 3.72, 3.91
    ]
group = rep(1:3, [4, 5, 3]);
exactkw(data, group)
# カイ二乗値 =  5.54872, 自由度 =  2, P 値 =  0.0623895
# 正確な P 値 =  0.0538961039
# 査察した分割表の個数は 27720 個

二重配列で与える場合 Array{Array{Float64,1},1}

x = [[3.42, 3.84, 3.96, 3.76], [3.17, 3.63, 3.47, 3.44, 3.39], [3.64, 3.72, 3.91]]
exactkw(x)
# カイ二乗値 =  5.54872, 自由度 =  2, P 値 =  0.0623895
# 正確な P 値 =  0.0538961039
# 査察した分割表の個数は 27720 個

頑張ってます。クリックお願いします。

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

Julia に翻訳--134　マン・ホイットニーの U 検定，exact test

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

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

マン・ホイットニーの U 検定(exact test)
http://aoki2.si.gunma-u.ac.jp/R/exact-mw.html

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

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

global は不要だった。

==========#

using Rmath, Printf

function exactmw(x, y = NaN; exact = true)
    function found()
        hh = abs(temp2 - sum(um[1, :] .* r))
        if hh >= stat_val || abs(hh - stat_val) <= EPSILON
            nprod = sum(perm_table[rt .+ 1]) - sum(perm_table[um .+ 1])
            nntrue += exp(nprod - nntrue2 * log_expmax)
            while nntrue >= EXPMAX
                nntrue /= EXPMAX
                nntrue2 += 1
            end
        end
        ntables += 1
    end

    function search(x, y)
        if y == 1
            found()
        elseif x == 1
            t = um[1, 1] - um[y, 1]
            if t >= 0
                um[1, 1] = t
                search(nc, y - 1)
                um[1, 1] += um[y, 1]
            end
        else
            search(x - 1, y)
            while um[y, 1] != 0 && um[1, x] != 0
                um[y, x] += 1
                um[y, 1] -= 1
                um[1, x] -= 1
                search(x - 1, y)
            end
            um[y, 1] += um[y, x]
            um[1, x] += um[y, x]
            um[y, x] = 0
        end
    end

    function exacttest()
        denom2 = 0
        denom = perm_table[n + 1] - sum(perm_table[ct .+ 1])
        while denom > log_expmax
            denom -= log_expmax
            denom2 += 1
        end
        denom = exp(denom)
        um[:, 1] = rt
        um[1, :] = ct
        search(nc, nr)
        @printf("正確な P 値 = % .10g\n", nntrue / denom * EXPMAX ^ (nntrue2 - denom2))
        @printf("査察した分割表の個数は % s 個\n", ntables)
    end

    function asymptotic()
        z = stat_val / sqrt(n1n2 / (n * (n - 1)) * (n ^ 3 - n - sum(ct .^ 3 .- ct)) / 12)
        @printf("U = % g, Z = % g, P 値 = % g\n", n1n2 / 2 - stat_val, z, pnorm(z, false) * 2)
    end

    if ndims(x) == 2
        t = x
    else
        nx = length(x)
        ny = length(y)
        indexx, indexy, t = table(rep(1:2, [nx, ny]), vcat(x, y))
    end
    EPSILON = 1e-10
    EXPMAX = 1e100
    log_expmax = log(EXPMAX)
    nr, nc = size(t)
    nr == 2 || error("nrow(table) wasn't 2")
    um = zeros(Int, nr, nc)
    rt = sum(t, dims=2)
    ct = transpose(sum(t, dims=1))
    n1 = rt[1]
    n2 = rt[2]
    n1n2 = n1 * n2
    n = n1 + n2
    r = cumsum(vcat(0, ct[1:nc-1])) .+ (ct .+ 1) ./ 2
    temp = n1n2 + n1 * (n1 + 1) / 2
    temp2 = temp - n1n2 / 2
    stat_val = abs(temp2 - sum(t[1, :] .* r))
    asymptotic()
    if exact
        perm_table = cumsum(vcat(0, log.(1:n .+ 1)))
        ntables = nntrue = nntrue2 = 0
        exacttest()
    end
end

function rep(x, n::Array{Int64,1})
    length(x) == length(n) || error("length(x) wasn't length(n)")
    vcat([repeat([x[i]], n[i]) for i = 1:length(x)]...)
end

function table(x, y) # 二次元
    indicesx = sort(unique(x))
    indicesy = sort(unique(y))
    counts = zeros(Int, length(indicesx), length(indicesy))
    for (i, j) in zip(indexin(x, indicesx), indexin(y, indicesy))
        counts[i, j] += 1
    end
    return indicesx, indicesy, counts
end

x = [5 3 2 1
     2 3 1 2];
exactmw(x)
# U =  33.5, Z =  0.907299, P 値 =  0.364248
# 正確な P 値 =  0.3837421608
# 査察した分割表の個数は 91 個

2 つのデータベクトルが与えられる場合

x = [78, 64, 75, 45, 82];
y = [110, 70, 53, 53];
exactmw(x, y)
[1 0 1 0 1 1 1 0; 0 2 0 1 0 0 0 1]
# U =  9, Z =  0.245976, P 値 =  0.805701
# 正確な P 値 =  0.8492063492
# 査察した分割表の個数は 91 個

頑張ってます。クリックお願いします。

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

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

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

プロフィール

プロフィール画像

バックナンバー

カレンダー

2021年3月
日 月 火 水 木 金 土
1 2 3 4 5 6
7 8 9 10 11 12 13
14 15 16 17 18 19 20
21 22 23 24 25 26 27
28 29 30 31

カテゴリー

雨雲の動き

  • RSS2.0