Julia で統計解析 一応の完

Julia で統計解析

2022/8/14 にバージョンアップ

https://r-de-r.github.io/stats/Julia-stats1.html
https://r-de-r.github.io/stats/Julia-stats2.html
https://r-de-r.github.io/stats/Julia-stats3.html
https://r-de-r.github.io/stats/Julia-stats4.html
https://r-de-r.github.io/stats/Julia-stats5.html
https://r-de-r.github.io/stats/Julia-stats6.html
https://r-de-r.github.io/stats/Julia-stats7.html
https://r-de-r.github.io/stats/Julia-stats8.html
https://r-de-r.github.io/stats/Julia-stats9.html

 

第 1 章　Julia の実行環境

1.　必要なファイルをダウンロードする
2.　Julia のインストール
3.　Julia を起動し終了する
4.　作業ディレクトリを変える
5.　Julia の環境設定をする
5.1.　startup.jl ファイルがあるかを確認
5.2.　config がなかった場合
5.3.　startup.jl がある場合
5.4.　変更内容の確認
6.　オンラインヘルプを使う
7.　パッケージを利用する
8.　エディタを使う
8.1.　Atom
8.2.　Jupyter lab
9.　結果を保存する

第 2 章　データの入出力

1.　データをデータフレームへ読み込む
1.1.　Excel のワークシートファイル
1.2.　CSV ファイル
1.2.1.　典型的な CSV ファイル
1.2.2.　一般の CSV ファイル
1.2.3.　CSV.read の引数
1.3. インターネット上の CSV ファイル
1.4.　モニター上に表示された表をデータフレームに読み込む
1.5.　クリップボードにコピーした内容をデータフレームにする
1.6.　文字列行から入力する
1.7.　エンコーディングの違うファイルを読む
1.8.　デリミタで区切られていない固定書式ファイルを読む
2.　データフレームを CSV ファイルに書き出す
2.1.　CSV.write の引数

第 3 章　データフレームの取り扱い

1.　データを使用するための準備
1.1.　既存のデータを使用する
1.2.　自前のデータ
2.　 データフレームの概要
2.1.　データフレームの大きさ
2.2.　データフレームの変数名
2.3.　データフレームの表示
3.　データフレームのコピーは copy()  で
4.　空データフレーム
5.　データフレームの列の参照
6.　データフレームの行の参照
7.　データフレームの列名の変更
8.　列の抽出
9.　指定された列を削除する
10.　行の抽出
11.　行の削除
12.　重複を除き，ユニークな行のみを含むデータフレームを作る
13.　欠損値を含む行か含まない行か
14.　欠損値を含まない行を抽出する
15.　データフレームの列方向連結
16.　データフレームの行方向連結
16.1.　cols = :setequal, cols=:orderequal
16.2.　cols = :intersect
16.3.　cols = :subset
16.4.　cols = :union
17.　データフレームの最終行の次に 1 行追加する
18.　行を繰り返してデータフレームを作る
19.　データフレームに列を挿入する
20.　データフレームの要素から新しいデータフレームを作る
21.　データフレームの各列に関数を施す
22.　ソート（並べ替え）
22.1.　ソートされているかをチェックする
22.2.　ソートの順番を決める
22.3.　並べ替えベクトルを返す
23.　行の包含
23.1.　df1 と df2 のすべての組み合わせを作る crossjoin
23.2.　df1 の行のうち，df2 に含まれない行を抽出する antijoin
23.3.　df1 に，df2 をマージする innerjoin
23.4.　df1 に，df2 をマージする leftjoin
23.5.　df1 に，df2 をマージする outerjoin
23.6.　df2 に，df1 をマージする rightjoin
23.7.　両方に存在する項目のみでマージする semijoin
24.　欠損値のリストワイズ除去
25.　ロングフォーマット（縦長データフレーム）に変換する
26.　ワイドフォーマット（横長データフレーム）に変換する
27.　データフレームをグループ変数に基づいて分割する
27.1. グループ分けされたデータフレームを抽出する
27.2.　グループ化されたデータフレームに関する情報
27.3.　親データフレームを返す
28.　列ごとに関数を適用する
29.　行ごとに関数を適用する
30.　クエリーによる操作例
30.1.　要素に関数を適用する　@map コマンド
30.2.　条件を満たす行を抽出　@filter コマンド
30.3.　データフレームのグループ化　@groupby コマンド
30.4.　データのソート @orderby，@orderby_descenidng，@thenby，@thenby_descending コマンド
30.5.　データフレームのマージ　@groupjoin コマンド
30.6.　データフレームの連結　@join コマンド
30.7.　キーと値のペアを展開　@mapmany コマンド
30.8.　要素を取り出す　@take コマンド
30.9.　要素を捨てる　@drop コマンド
30.10.　重複データを除く　@unique コマンド
30.11.　列の選択　@select コマンド
30.12.　列名の変更　@rename コマンド
30.13.　変数変換　@mutate コマンド
30.14.　欠損値行を除く　@dropna コマンド
30.15.　欠損値を含まないデータフレーム　@dissallowna コマンド
30.16.　欠損値の置き換え　@replacena コマンド
31.　データフレームを二次元配列に変換する

第 4 章　一変量統計

1.　一変量統計
1.1.　一変数の場合
1.1.1.　基礎統計量
1.1.2.　パーセンタイル値
1.1.3.　度数分布
1.2.　複数の変数の場合
1.2.1.　eachcol() を使う
1.2.2.　describe() を使う
1.2.3.　combine(), select()/select!(), transform()/transform!() を使う
1.3.　グループごとの記述統計量
1.3.1.　describe() を使う
1.3.2.　変数ごとに統計量を一覧表示
2.　二変量統計
2.1.　二重クロス集計表
2.2.　相関係数，共分散
2.2.1.　欠損値を含まない場合
2.2.1.1.　それぞれの関数を使う
2.2.1.2.　combine() を使う
2.2.2.　欠損値を含む場合
2.2.2.1.　それぞれの関数を使う
2.2.2.2.　combine() を使う
2.3.　相関係数行列，分散・共分散行列
2.3.1.　欠損値を含まない場合
2.3.2.　欠損値を含む場合
3.　多変量統計
3.1. 多重クロス集計表

第 5 章　離散データの可視化

1.　例示に使用するデータセット
1.1.　カテゴリーデータ
2.　棒グラフ
2.1.　一標本の場合
2.2.　二標本以上の場合
2.2.1.　横に並べる棒グラフ
2.2.2.　積み上げ棒グラフ
2.3.　複数のグラフを行列状にまとめて表示する方法
3.　帯グラフ
4.　モザイクプロット
5.　バルーンプロット

第 6 章　数値データの可視化

1.　ヒストグラム
1.1.　一標本の場合
1.2.　複数標本の場合
2.　ボックスプロット（箱ひげ図）
3.　バイオリンプロット
4.　ドットプロット
5.　カーネル密度推定の描画
6.　Q-Q プロット
7.　散布図
8.　カーネル密度推定
9.　散布図行列
10.　作図レシピ
10.1.　二軸グラフ
10.2.　Andrews plot
10.3.　星座グラフ
10.4.　サンフラワープロット
10.5.　レーダーチャート

第 7 章　検定と推定

1.　検定関数関数の呼び出し方
2.　検定関数関数により得られる結果の利用法
3.　分布の検定
3.1.　観察度数が一様かどうかの検定
3.1.1.　ピアソンの \\chi^2 検定
3.1.2.　対数尤度比検定（G^2 検定）
3.2.　観察度数が理論比に從うかどうかの検定
3.2.1.　ピアソンの \\chi^2 検定
3.2.2.　対数尤度比検定（G^2 検定）
4.　独立性の検定
4.1.　ピアソンの \\chi^2 検定
4.2.　対数尤度比検定（G^2 検定）
5.　パワーダイバージェンス検定
6.　フィッシャーの正確検定
7.　二項検定
8.　t 検定
8.1.　一標本の検定（母平均の検定）
8.2.　等分散の場合の t 検定
8.3.　等分散でない場合 Welch の方法
9.　マン・ホイットニーの U 検定
10.　符号検定
11.　ウィルコクソンの符号付き順位和検定
12.　相関係数の検定
13.　対応のない k 標本（独立 k 標本）
13.1.　一元元配置分散分析
13.2.　一元元配置分散分析（ウェルチの方法）
13.3.　クラスカル・ウォリス検定
14.　コルモゴロフ・スミルノフ検定
14.1.　1 標本の分布の検定
14.2.　2 標本の分布の差の検定
15.　アンダーソン・ダーリング検定
15.1.　1 標本の場合
15.2.　k 標本の場合
16.　ワルド・ウォルフォビッツ連検定
17.　並べ替え検定（無作為検定）

第 8 章　多変量解析

1. 回帰分析
1.1.　線形最小二乗回帰 Linear Least Square
1.1.1.　単回帰分析
1.1.2.　重回帰分析
1.2.　リッジ回帰 Ridge Regression
1.2.1.　説明変数が 1 個の場合
1.2.2.　説明変数が 2 個以上の場合
1.3.　GLM パッケージによる回帰分析
1.3.1.　重回帰分析 Ordinary Least Squares Regression
1.3.2.　プロビット回帰 Probit Regression
1.3.3.　ロジット回帰 Logit Regression
1.4.　多項式回帰
1.4.1.　重回帰プログラムを用いる
1.4.2.　多項式回帰 Polynomials パッケージ
1.5.　指数回帰
1.5.1.　説明変数が 1 個の場合
1.5.2.　説明変数が 2 個以上の場合
1.6.　累乗回帰
1.6.1.　説明変数が 1 個の場合
1.6.2.　説明変数が 2 個以上の場合
1.7.　非線形回帰
2.　判別分析
2.1.　二群判別分析
2.2.　多重判別分析
3.　主成分分析
4.　因子分析
5.　古典的多次元尺度解析
6.　クラスター分析
6.1.　K-means 法による非階層的クラスター分析
6.2. 階層的クラスター分析
7.　カテゴリー変数の取り扱い方
7.1.　重回帰分析の場合
7.2.　判別分析の場合

第 9 章　関数

 

 

「ナイチンゲールの鶏頭図よりは折れ線グラフ」を Julia で描く

ナイチンゲールの鶏頭図よりは折れ線グラフ...という R のプログラムも，あまり良くない気がして，Julia で書いて（描いて）みた。

R の example(Nightingale) の図

折れ線に重ねるデータポイントが，文字通り「間が抜けている」

Julia で描いた図。自画自賛。

using RDatasets, DataFrames
Nightingale = dataset("HistData", "Nightingale");

using Plots
function line(x, y, col, legend)
    plot!(x, y, color=col, label=legend)
    scatter!(x, y, markersize=3, color=col)
end

using Dates
function linechart(df)
    gr(grid=false, tick_direction=:out, titlefont=(11, "times"),
        guidefont=(10, "times"), tickfont=(8, "times"),
        markerstrokewidth=0, fontfamily=(8, "serif"),
        foreground_color_legend=nothing, label="")
    tmtick  = Date(1854, 4):Month(6):Date(1856, 3)
    tmticks = Dates.format.(tmtick, "u yyyy")
    plot(title="Causes of Mortality of the British Army in the East",
        xlabel="Month/Year", ylabel="Annual Death Rate",
        xticks=(tmtick, tmticks), xminorticks=6)
    line(df.Date, df.DiseaseRate, :blue, "Preventable disease")
    line(df.Date, df.WoundsRate, :red, "Wounds and injuries")
    line(df.Date, df.OtherRate, :black, "Other")
    plot!(df.Date[[1, 12, 12, 1]], [0, 0, 1100, 1100],
        seriestype=:shape, fillcolor="gray10", alpha=0.04, legend=:right)
    annotate!(df.Date[1], 1000, text("Before Sanitary Commission", "times", :left, 10))
    annotate!(df.Date[13], 1000, text("After Sanitary Commission", "times", :left, 10))
end

linechart(Nightingale)

ナイチンゲールの鶏頭図--再び

2019/12/10 にもナイチンゲールの鶏頭図について書いたが，それは R であった。

そのときには，HistData パッケージの Nightingale は知らなかった。

中澤先生が example(Nightingale) で描けると言及していたので，やってみた。

melt がないとか，reshape がないとか，ggplot2 がないとかうるさかった割に，ちょっと残念な出来栄えだった。

ので，Julia で描いたみた。ggplot の呪いがないので，スッキリ書けたと思う。プログラムは下の方に。

2 枚の図は描画範囲を一致させている。ちょっと，色がどぎつい。

using RDatasets, DataFrames
Nightingale = dataset("HistData", "Nightingale");

using Colors, Plots, Dates

function spidernet(df, begin_)
    θ = range(0, 360, length=1200)
    for deg in 250:250:1250
        r = sqrt(deg)
        plot!(r * cosd.(θ), r * sind.(θ), color=:gray70)
        annotate!(r * cosd(150), r * sind(150),
            text(string(deg), 6, :gray70, "times"))
    end
    for i in 1:12
        pos = sqrt(1200)
        x = pos*cosd(begin_[i] - 15)
        y = pos*sind(begin_[i] -15)
        plot!([0, x], [0, y], color=:gray70)
        annotate!(x, y,
            text(Dates.format.(df[i, 1], "yyyy/mm"), 8,
                :gray70, "times", 4<=i<=9 ? :left : :right))
    end
end

function legend_(df, color)
    for (y, col, leg) in zip(-14:-3:-20, color,
                             ["■ Disease", "■ Other", "■ Wounds"])
        annotate!(32, y, text(leg, 8, :left, col, "times"))
    end
end

function nightingale(df)
    begin_ = 180:-30:-210
    # color = [colorant"#b1bbc1", colorant"#e6b4b7", colorant"#6f625e"];
    color = [colorant"#f6716d", colorant"#6298ff", colorant"#1fbf38", ];
    gr(grid=false, showaxis=false, ticks=false, aspect_ratio=1, lw=0.2,
        titlefont = (10, "times"), label="")
    p = plot(title="Diagram of the Causes of Mortality in the Army in the East")
    spidernet(df, begin_)
    legend_(df, color)
    for i in 1:12
        deg = range(begin_[i], begin_[i+1], length = 100)
        for j in sortperm(Vector(df[i, 8:10]), rev = true)
            r = sqrt(df[i, j+7])
            x = vcat(0, r * cosd.(deg))
            y = vcat(0, r * sind.(deg))
            plot!(x, y, seriestype=:shape, fillcolor = color[j], label="")
        end
    end
    p
end

nightingale(Nightingale[1:12, :])
nightingale(Nightingale[13:24, :])

 

Julia: sum(), var() に対応する describe()

"""
DataFrames　の describe() は，sum(), var() を計算しない。

sum(), var() に対応する関数 describe2() を定義しておく。

指定できるのは :sum, :mean, :var, :std, :min, :q25, :median, :q75, :max, :nmissing である。

すべての統計量を計算するときは，統計量の指定をしないようにするか :all を指定すればよい。

なお，DataFrames.describe() と競合してもかまわないならば，関数名を describe にしてもよい。
この場合，DataFrames　の describe() を使うときは，DataFrames.describe() とする。

"""

using DataFrames, Statistics
function describe(df, arg...=:all)
    q25(x) = quantile(x, 0.25)
    q75(x) = quantile(x, 0.75)
    nmissing(x) = nrow(df) - length(ismissing.(x))
    arg0 = [:sum, :mean, :var, :std, :min, :q25, :median, :q75, :max, :nmissing]
    func = [sum, mean, var, std, minimum, q25, median, q75, maximum, nmissing]
    funcname = ["sum", "mean", "var", "std", "minimum", "q25", "median", "q75", "maximum", "nmissing"]
    arg[1] == :all && (arg = arg0)
    res = DataFrame(variable = Symbol.(names(df)))
    n = ncol(df)
    m = length(arg)
    body = Array{Float64, 2}(undef, n, m)
    loc = indexin(arg, arg0)
    for i in 1:n
        x = skipmissing(df[!, i])
        body[i, :] = [func[j](x) for j in loc]
    end
    bodydf = DataFrame(body, funcname[loc])
    pos = indexin(["nmissing"], funcname[loc])[1]
    isnothing(pos) || (bodydf[!, pos] = Int.(bodydf[:, pos]))
    hcat(res, bodydf)
end

"""
使用例

julia> using RDatasets, Statistics

julia> iris = dataset("datasets", "iris");

julia> describe(iris[:, 1:4])
4×11 DataFrame
 Row │ variable     sum      mean     var       std       minimum  q25      median   q75      maximum  nmissing
     │ Symbol       Float64  Float64  Float64   Float64   Float64  Float64  Float64  Float64  Float64  Int64
─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ SepalLength    876.5  5.84333  0.685694  0.828066      4.3      5.1     5.8       6.4      7.9         0
   2 │ SepalWidth     458.6  3.05733  0.189979  0.435866      2.0      2.8     3.0       3.3      4.4         0
   3 │ PetalLength    563.7  3.758    3.11628   1.7653        1.0      1.6     4.35      5.1      6.9         0
   4 │ PetalWidth     179.9  1.19933  0.581006  0.762238      0.1      0.3     1.3       1.8      2.5         0

julia> DataFrames.describe(iris[:, 1:4])
4×7 DataFrame
 Row │ variable     mean     min      median   max      nmissing  eltype
     │ Symbol       Float64  Float64  Float64  Float64  Int64     DataType
─────┼─────────────────────────────────────────────────────────────────────
   1 │ SepalLength  5.84333      4.3     5.8       7.9         0  Float64
   2 │ SepalWidth   3.05733      2.0     3.0       4.4         0  Float64
   3 │ PetalLength  3.758        1.0     4.35      6.9         0  Float64
   4 │ PetalWidth   1.19933      0.1     1.3       2.5         0  Float64
"""

Julia で水平方向のデンドログラムを描く

"""

Julia の Clustering.hclust() は StatsPlots パッケージで普通の（垂直方向の）デンドログラムを描く。機能を持たない。

以下のプログラムは，R の plot() と同じであるが水平方向のデンドログラムを描く。
"""

using Plots
function plot_hclust_horizontal(hc)
    function get(i, x)
        x[i] < 0 && return(ord[abs(x[i])])
        get(x[i], x)
    end
    n = length(hc.order)
    apy = collect(n:-1:1) .+ 0.5
    apx = zeros(n)
    ord = sortperm(hc.order)
    plot(yshowaxis=false, yticks=false, tick_direction=:out, grid=false,
         xlims=(-0.05, maximum(hc.height)), ylims=(0, apy[1]),
         xlabel="height", label="")
    for i in 1:n
        annotate!(0, apy[i], text(string(hc.labels[hc.order[i]]) * " ", 8, :right))
    end
    for i in 1:n-1
        c1 = get(i, hc.merge[:,1])
        c2 = get(i, hc.merge[:,2])
        plot!([apx[c1], hc.height[i], hc.height[i], apx[c2]],
              [apy[c1], apy[c1], apy[c2], apy[c2]], color=:black, label="")
        apx[c1] = apx[c2] = hc.height[i]
        apy[c1] = apy[c2] = (apy[c1] + apy[c2]) / 2
    end
    plot!()
end

"""
使用法

julia> using RDatasets

julia> iris = RDatasets.dataset("datasets", "iris");

julia> function distancematrix(X)
           nr,  nc = size(X)
           d = zeros(nr, nr)
           for i in 1:nr-1
               for j in i+1:nr
                   d[i, j] = d[j, i] = sqrt(sum((X[i, :] .- X[j, :]).^2))
               end
           end
           d
       end;

julia> # using Random; Random.seed!(123) # 毎回同じ結果にするためには乱数の種を設定する
       using Clustering

julia> x = Matrix(iris[1:20, 1:4]);

julia> D = distancematrix(x);

julia> hc = hclust(D, linkage=:ward);

julia> plot_hclust_horizontal(hc)

"""

Julia で統計解析--その8 多変量解析

Julia で統計解析--その8　多変量解析

これらの文書群は github で管理することとした

最新バージョン 2022-02-22 07:25

以下を参照のこと

https://r-de-r.github.io/stats/Julia-stats8.html

1. 回帰分析
1.1.　線形最小二乗回帰 Linear Least Square
1.1.1.　単回帰分析
1.1.2.　重回帰分析
1.2.　リッジ回帰 Ridge Regression
1.2.1.　説明変数が 1 個の場合
1.2.2.　説明変数が 2 個以上の場合
1.3.　GLM パッケージによる回帰分析
1.3.1.　重回帰分析 Ordinary Least Squares Regression
1.3.2.　プロビット回帰 Probit Regression
1.3.3.　ロジット回帰 Logit Regression
1.4.　多項式回帰
1.4.1.　重回帰プログラムを用いる
1.4.2.　多項式回帰 Polynomials パッケージ
1.5.　指数回帰
1.5.1.　説明変数が 1 個の場合
1.5.2.　説明変数が 2 個以上の場合
1.6.　累乗回帰
1.6.1.　説明変数が 1 個の場合
1.6.2.　説明変数が 2 個以上の場合
1.7.　非線形回帰
2.　判別分析
2.1.　二群判別分析
2.2.　多重判別分析
3.　主成分分析
4.　因子分析
5.　古典的多次元尺度解析
6.　クラスター分析
6.1.　K-means 法による非階層的クラスター分析
6.2. 階層的クラスター分析
7.　カテゴリー変数の取り扱い方
7.1.　重回帰分析
7.2.　判別分析の場合

1 行できれいな散布図描けるのに，なんで ggplot2?

using Plots
z = randn(1000, 100);
x = sum(z[:, 1:70], dims=2);
y = sum(z[:, 30:100], dims=2);
scatter(x, y,
    grid=false, # なくてもいい
    tick_direction=:out, # なくてもいい
    label="",  # なくてもいい
    )

1 行できれいなヒストグラム描けるのに，なんで ggplot2?

using Plots
x = randn(10000)
histogram(x,
    grid=false, # なくてもいい
    tick_direction=:out, # なくてもいい
    label="",  # なくてもいい
    )

Julia で統計解析--その3 データフレームの取り扱い

これらの文書群は github で管理することとした

最新バージョン 2022-02-12 22:38

以下を参照のこと

https://r-de-r.github.io/stats/Julia-stats3.html

1.　データを使用するための準備
 1.1.　既存のデータを使用する
 1.2.　自前のデータ
2.　 データフレームの概要
 2.1.　データフレームの大きさ
 2.2.　データフレームの変数名
 2.3.　データフレームの表示
3.　データフレームのコピーは copy()  で
4.　空データフレーム
5.　データフレームの列の参照
6.　データフレームの行の参照
7.　データフレームの列名の変更
8.　列の抽出
9.　指定された列を削除する
10.　行の抽出
11.　行の削除
12.　重複を除き，ユニークな行のみを含むデータフレームを作る
13.　欠損値を含む行か含まない行か
14.　欠損値を含まない行を抽出する"
15.　データフレームの列方向連結
16.　データフレームの行方向連結
 16.1.　cols = :setequal, cols=:orderequal
 16.2.　cols = :intersect
 16.3.　cols = :subset
 16.4.　cols = :union
17.　データフレームの最終行の次に 1 行追加する
18.　行を繰り返してデータフレームを作る
19.　データフレームに列を挿入する
20.　データフレームの要素から新しいデータフレームを作る
21.　データフレームの各列に関数を施す
22.　ソート（並べ替え）
 22.1.　ソートされているかをチェックする
 22.2.　ソートの順番を決める
 22.3.　並べ替えベクトルを返す
23.　行の包含
 23.1.　df1 と df2 のすべての組み合わせを作る crossjoin
 23.2.　df1 の行のうち，df2 に含まれない行を抽出する antijoin
 23.3.　df1 に，df2 をマージする innerjoin
 23.4.　df1 に，df2 をマージする leftjoin
 23.5.　df1 に，df2 をマージする outerjoin
 23.6.　df2 に，df1 をマージする rightjoin
 23.7.　両方に存在する項目のみでマージする semijoin
24.　欠損値のリストワイズ除去
25.　ロングフォーマット（縦長データフレーム）に変換する
26.　ワイドフォーマット（横長データフレーム）に変換する
27.　データフレームをグループ変数に基づいて分割する
 27.1. グループ分けされたデータフレームを抽出する
 27.2.　グループ化されたデータフレームに関する情報
 27.3.　親データフレームを返す
28.　列ごとに関数を適用する
29.　行ごとに関数を適用する
30.　クエリーによる操作例
 30.1.　要素に関数を適用する　@map コマンド
 30.2.　条件を満たす行を抽出　@filter コマンド
 30.3.　データフレームのグループ化　@groupby コマンド
 30.4.　データのソート @orderby，@orderby_descenidng，@thenby，@thenby_descending コマンド
 30.5.　データフレームのマージ　@groupjoin コマンド
 30.6.　データフレームの連結　@join コマンド
 30.7.　キーと値のペアを展開　@mapmany コマンド
 30.8.　要素を取り出す　@take コマンド
 30.9.　要素を捨てる　@drop コマンド
 30.10.　重複データを除く　@unique コマンド
 30.11.　列の選択　@select コマンド
 30.12.　列名の変更　@rename コマンド
 30.13.　変数変換　@mutate コマンド
 30.14.　欠損値行を除く　@dropna コマンド
 30.15.　欠損値を含まないデータフレーム　@dissallowna コマンド
 30.16.　欠損値の置き換え　@replacena コマンド
31.　データフレームを二次元配列に変換する

 

Julia で統計解析--その7 検定と推定

これらの文書群は github で管理することとした

最新バージョン 2022-02-12 22:23

以下を参照のこと

https://r-de-r.github.io/stats/Julia-stats7.html

1. HypothesisTests に含まれる検定関数の使用法
 1.1.　検定関数関数の呼び出し方
 1.2.　検定関数関数により得られる結果の利用法
2.　検定と推定
 2.1.　分布の検定
  2.1.1.　観察度数が一様かどうかの検定
   2.1.1.1.　ピアソンのχ二乗検定
   2.1.1.2.　対数尤度比検定（G2 検定）
  2.1.2.　観察度数が理論比に從うかどうかの検定
   2.1.2.1.　ピアソンのχ二乗検定
   2.1.2.2.　対数尤度比検定（G2 検定）
 2.2.　独立性の検定
  2.2.1.　ピアソンのχ二乗検定
  2.2.2.　対数尤度比検定（G2 検定）
 2.3.　パワーダイバージェンス検定
 2.4.　フィッシャーの正確検定
 2.5.　二項検定
 2.6.　t 検定
  2.6.1.　一標本の検定（母平均の検定）
  2.6.2.　等分散の場合の t 検定
  2.6.3.　等分散でない場合 Welch の方法
 2.7.　マン・ホイットニーの U 検定
 2.8.　符号検定
 2.9.　ウィルコクソンの符号付き順位和検定
 2.10.　相関係数の検定
 2.11.　対応のない k 標本（独立 k 標本）
  2.11.1.　一元元配置分散分析
  2.11.2.　一元元配置分散分析（ウェルチの方法）
  2.11.3.　クラスカル・ウォリス検定
 2.12.　コルモゴロフ・スミルノフ検定
  2.12.1.　1 標本の分布の検定
  2.12.2.　2 標本の分布の差の検定
 2.13.　アンダーソン・ダーリング検定
  2.13.1.　1 標本の場合
  2.13.2.　k 標本の場合
 2.14.　ワルド・ウォルフォビッツ連検定
 2.15.　並べ替え検定（無作為検定）

 

Julia が 1.7.2 にアップデートされた

うれしいなぁ

https://julialang.org/downloads/#current_stable_release

Julia: RCall が Julia 1.7.2 / macOS Monterey / M1 チップで 動くようになった

しばらく前からというのも，Julia 1.7.1 になったり，mac OS Monterey になったり，マシンが M1 チップ搭載の Mac mini になったりで，わたしの環境で RCall が動かなくなっていた。

今日，たまたまなんの気無しにやってみたら，動いた。ばんざい！！

julia> using RCall

julia> R"""
       x <- matrix(c(3, 5, 6, 8), byrow=TRUE, nc=2)
       fisher.test(x)
       """
RObject{VecSxp}

 Fisher's Exact Test for Count Data

data:  x
p-value = 1
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.08839317 6.37801297
sample estimates:
odds ratio 
 0.8081216 

Julia が 1.7.2 になったせいかな?

Current stable release: v1.7.2 (Feb 6, 2022)

Julia: 一般項がn、n²、n³、...、n^10、n^50の数列の和

SymPy で n^10 までと，n^50 をやってみました

julia> summation(k^50, (k, 1, n)) |> factor |> print

n*(n + 1)*(2*n + 1)*(429*n^48 + 10296*n^47 + 75504*n^46 - 118404*n^45 - 3433716*n^44 + 5209776*n^43 + 177855964*n^42 - 269388834*n^41 - 8790556346*n^40 + 13320528936*n^39 + 399690120544*n^38 - 606195445284*n^37 - 16528902129916*n^36 + 25096450917516*n^35 + 617619215226134*n^34 - 938977048297959*n^33 - 20729072134582867*n^32 + 31563096726023280*n^31 + 620996836434449376*n^30 - 947276803014685704*n^29 - 16489621590673441320*n^28 + 25208070787517504832*n^27 + 385039595475766105512*n^26 - 590163428607407910684*n^25 - 7834990698171041513900*n^24 + 12047567761560266226192*n^23 + 137484855642671592627616*n^22 - 212251067344787522054520*n^21 - 2054997867949212429383888*n^20 + 3188622335596212405103092*n^19 + 25783333201061356027601406*n^18 - 40269310969390140243953655*n^17 - 266742505880274126845610963*n^16 + 420248414305106260390393272*n^15 + 2225375820438985887388497936*n^14 - 3548187937811031961277943540*n^13 - 14547150581890312152150257764*n^12 + 23594819841740984208864358416*n^11 + 71655262530103266642836795708*n^10 - 119280303716025392068687372770*n^9 - 251166419012752676999004054970*n^8 + 436389780377141711532849768840*n^7 + 569535901810418645522117212000*n^6 - 1072498742904198824049600702420*n^5 - 679511742411895276079901986300*n^4 + 1555516985069942326144653330660*n^3 + 137940192459690838148194087930*n^2 - 984668781224507420294617797225*n + 328222927074835806764872599075)/43758

Julia: sin7.5°、cos7.5°、tan7.5°はどんな数？

思ったより，簡単で，びっくりした。

Julia: HypothesisTests.BinomialTest() の p 値に難あり

HypothesisTests の二項検定の p 値が駄目だ。正確に言うと，母比率 ≒ 0.5 のときの両側検定の場合に限るのだが。

母比率 ≒ 0.5のときは，確率分布が歪んでいる。

例えば，x = 18, n = 24, p = 0.68 だと，以下のようになる。

julia> using Plots

julia> bar(pdf(obj), xticks=(1:25, 0:24), xlimit=(7, 26), grid=false, tick_direction=:out, label="")

julia> using DataFrames

julia> obj = Binomial(24, 0.68)
Binomial{Float64}(n=24, p=0.68)

julia> p = [pdf(obj, x) for x in 0:24];

julia> cum1 = cumsum(p);

julia> cum2 = reverse(cumsum(reverse(p)));

julia> df = DataFrame(x = 0:24, p = p, cum1 = cum1, cum2 = cum2)
25×4 DataFrame
 Row │ x      p            cum1         cum2       
     │ Int64  Float64      Float64      Float64    
─────┼─────────────────────────────────────────────
   1 │     0  1.32923e-12  1.32923e-12  1.0
   2 │     1  6.77906e-11  6.91199e-11  1.0
   3 │     2  1.65663e-9   1.72575e-9   1.0
   4 │     3  2.58159e-8   2.75416e-8   1.0
   5 │     4  2.88008e-7   3.1555e-7    1.0
   6 │     5  2.44807e-6   2.76362e-6   1.0
   7 │     6  1.64735e-5   1.92371e-5   0.999997
   8 │     7  9.00158e-5   0.000109253  0.999981
   9 │     8  0.000406477  0.00051573   0.999891
  10 │     9  0.00153558   0.00205131   0.999484
  11 │    10  0.00489467   0.00694598   0.997949
  12 │    11  0.0132378    0.0201838    0.993054
  13 │    12  0.0304746    0.0506585    0.979816
  14 │    13  0.0597772    0.110436     0.949342
  15 │    14  0.0998065    0.210242     0.889564
  16 │    15  0.141393     0.351635     0.789758
  17 │    16  0.169008     0.520643     0.648365
  18 │    17  0.169008     0.689651     0.479357
  19 │    18  0.139667     0.829318     0.310349
  20 │    19  0.0937236    0.923041     0.170682
  21 │    20  0.0497907    0.972832     0.0769586
  22 │    21  0.0201534    0.992985     0.0271679
  23 │    22  0.0058399    0.998825     0.00701455
  24 │    23  0.00107911   0.999904     0.00117466
  25 │    24  9.55463e-5   1.0          9.55463e-5

x = 18 が平均値 n * p = 24 * 0.68 = 16.32 より大きいので，18〜24 の確率の合計を求めて「2倍してしまう」これは大きな間違い（大したことではないという人もいるんだが）。

正しくは，x = 18 と平均値を挟んで反対側にあるx = 18 のときの確率より大きいものの最小値から1引いた14 から 0 までの確率の和，青で色付した 0.210242 と，x = 18~24 の確率の和，赤で色付した 0.310349 の和，0.210242 + 0.310349 = 0.520591 が求める p 値である。

ということで，以下のプログラムになってしまった。

最初っから書いてももっと短くなるかもしれないが，どういうわけか母比率の信頼区間は正しく計算されているようなので，その部分は残した(confint をそのまま使って）。

使い方と結果の出力を，R に似せておいた。

binom_test(18, 24, p = 0.68, conf_level=0.98, alternative="two.sided")

Exact binomial test

number of successes = 18, number of trials = 24, p-value = 0.52059
alternative hypothesis: true probability of success is not equal to 0.68
98 percent confidence interval: [0.4951508, 0.9200126]
sample estimates: nprobability of success = 0.75

#####

using HypothesisTests
using Distributions

roundn(x, digits) = round(x, digits=digits)
formatp(p) = p < 0.00001 ? "< 0.00001" : string(roundn(p, 5))

function alt_type(alternative)
    res = indexin([alternative], ["two.sided", "less", "greater"])[1]
    if isnothing(res)
        println("alternative is invalid, and changed to 'two.sided'")
        res = 1
    end
    res
end

function binom_test(x, n; p=0.5, conf_level=0.95, alternative="two.sided")
    function pvalue(x, n, p, tail)
        obj = Binomial(n, p)
        if tail == :left
            cdf(obj, x)
        elseif tail == :right
            ccdf(obj, x - 1)
        else
            if p == 0
                Float64(x == 0)
            elseif p == 1
                Float64(x == n)
            else
                limit = 1.0000001pdf(obj, x)
                m = n * p
                if x == m
                    1
                elseif x < m
                    y = sum(pdf.(obj, ceil(Int, m):n) .<= limit)
                    cdf(obj, x) + ccdf(obj, n - y)
                else
                    y = sum(pdf.(obj, 0:floor(Int, m)) .<= limit)
                    cdf(obj, y - 1) + ccdf(obj, x - 1)
                end
            end
        end
    end
    ###
    type = alt_type(alternative)
    word = ["not equal to", "less than", "greater than"][type]
    tail = [:both, :left, :right][type]
    result = BinomialTest(x, n, p);
    p_value = pvalue(x, n, p, tail)
    lo, hi = confint(result, level=conf_level, tail=tail)
    println("\tExact binomial test\n\ndata:  x and n")
    println("number of successes = $x, number of trials = $n, p-value = $(formatp(p_value))")
    println("alternative hypothesis: true probability of success is $word $p")
    println("$(round(Int, 100conf_level)) percent confidence interval: [$(roundn(lo, 7)), $(roundn(hi, 7))]")
    println("sample estimates: nprobability of success = $(roundn(x/n, 7))\n\n")
end