Help us understand the problem. What is going on with this article?

Juliaで遊んでみた - 確率分布・仮説検定と可視化編

最近流行りつつある科学計算処理向けの高水準言語「Julia」の基本操作をまとめていきます。
1.準備編
2.確率分布・仮説検定と可視化編
3.データフレーム操作編

chapter.2 確率分布

確率分布のパッケージがあったので、ほぼそのパッケージの紹介のような感じです。
R とはまた使用感が異なりますが、個人的にはこっちの方が直感的で好きです。

パッケージ:Distributions.jl
詳細:https://juliastats.org/Distributions.jl/stable/

1.パッケージ install

julia> import Pkg; Pkg.add("Distributions")
julia> using Distributions

この後、ランダムにデータを生成する場面がありますが、seed を固定したい場合は以下

julia> using Random
julia> Random.seed!([シード値])

2.確率分布オブジェクト

とりあえず名称そのま関数を実行してみるとわかる通り、デフォルトのパラメーターが設定されている。もちろん指定することもできる。

julia> Bernoulli()
Bernoulli{Float64}(p=0.5)
julia> Binomial()
Binomial{Float64}(n=1, p=0.5)
julia> Poisson()
Bernoulli{Float64}(p=0.5)

julia> Normal()
Normal{Float64}(μ=0.0, σ=1.0)
julia> Exponential()
Exponential{Float64}(θ=1.0)
julia> Gamma()
Gamma{Float64}(α=1.0, θ=1.0)

R で言う所の rnorm,pnrom,dnorm,qnormをやってみる。

julia> norm_dist = Normal(50.0, 10.0) # 確率分布オブジェクトを生成
Normal{Float64}(μ=50.0, σ=10.0)

julia> rand(norm_dist, 5) # 第二引数に size
5-element Array{Float64,1}:
 31.32159859451268 
 76.52740439407731 
 43.67613857974404 
 49.230729489493726
 38.79220085610626

julia> rand(norm_dist, 5) # 第三引数で次元数
5×2 Array{Float64,2}:
 53.1482  33.9218
 40.4533  62.3337
 44.2832  45.7319
 48.4506  53.5824
 46.6584  39.3576

julia> pdf(norm_dist, 70)
0.005399096651318806

julia> cdf(norm_dist, 70)
0.9772498680518208

julia> quantile(norm_dist, 0.5)
50.0

パラメータを取り出す

julia> println("μ : ", norm_dist.μ)
μ : 50.0

julia> println("σ : ", norm_dist.σ)
σ : 10.0

配列の場合

julia> norm_dist = Normal() # とりあえずデータを生成
Normal{Float64}(μ=0.0, σ=1.0)

julia> array = rand(norm_dist, 5)
5-element Array{Float64,1}:
 -0.615182929965525  
  0.3018546554695042 
  0.43501627608696825
 -1.5087403568012359 
 -1.5281565440128444 

関数名称の直後に . を入れるだけ

julia> println(pdf.(norm_dist, array))
[0.33016473508059785, 0.3811750160113275, 0.36292534993233966, 0.1278257506789366, 0.12411214279395513]

julia> println(cdf.(norm_dist, array))
[0.26921695979105614, 0.6186185680364118, 0.6682246933064511, 0.06568257397844732, 0.06323683764142299]

julia> println(quantile.(norm_dist, 0.025 : 0.025 : 0.1))
[-1.9599639845400592, -1.6448536269514724, -1.439531470938456, -1.2815515655446004]

最尤推定

julia> norm_fitted = fit(Normal, array) # ここでfit
Normal{Float64}(μ=-0.5830417798446266, σ=0.8450652913045483)

確認してみる

julia> println(norm_fitted)
Normal{Float64}(μ=-0.5830417798446266, σ=0.8450652913045483)

パラメータを取り出す場合は先ほど同様...

julia> println(norm_fitted.μ)
-0.5830417798446266

julia> println(norm_fitted.σ)
0.8450652913045483

chapter.3 仮説検定

こちらも仮説検定パッケージがあったのでほぼその紹介になります。
パッケージ:HypothesisTests
詳細:https://hypothesistestsjl.readthedocs.io/en/latest/index.html

1.パッケージ install

julia> import Pkg; Pkg.add("HypothesisTests")
julia> using HypothesisTests, Distributions

2.とりあえずT検定

一標本

julia> x = rand(Normal(1,1), 10) # 適当にデータ生成

julia> ttest_object = OneSampleTTest(x, 0) # 第一引数にデータ、第二引数に帰無仮説
One sample t-test
-----------------
Population details:
    parameter of interest:   Mean
    value under h_0:         0
    point estimate:          1.332171049138298
    95% confidence interval: (0.7003, 1.9641)

Test summary:
    outcome with 95% confidence: reject h_0
    two-sided p-value:           0.0010

Details:
    number of observations:   10
    t-statistic:              4.7689293088672375
    degrees of freedom:       9
    empirical standard error: 0.27934384488805275

デフォルトでは両側なので、片側検定はpvalue関数を使う。tail で どっち側かを指定する

julia> println("right-side p-value : ", pvalue(ttest_object, tail=:right))
right-side p-value : 0.0005084438096162109

一応 Paired T

julia> before = rand(Normal(0,3), 10)
julia> after  = rand(Normal(1,1), 10)
julia> diff = after - before
julia> OneSampleTTest(diff, 0)
One sample t-test
-----------------
Population details:
    parameter of interest:   Mean
    value under h_0:         0
    point estimate:          0.661240160838284
    95% confidence interval: (-1.4651, 2.7876)

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.4996

Details:
    number of observations:   10
    t-statistic:              0.70347237950898
    degrees of freedom:       9
    empirical standard error: 0.9399660599323403

等分散を仮定したT検定

julia> x1 = rand(Normal(0,1), 10)
julia> x2 = rand(Normal(1,1), 20)
julia> println(EqualVarianceTTest(x1, x2))
Two sample t-test (equal variance)
----------------------------------
Population details:
    parameter of interest:   Mean difference
    value under h_0:         0
    point estimate:          -1.3055004613593664
    95% confidence interval: (-1.9006, -0.7104)

Test summary:
    outcome with 95% confidence: reject h_0
    two-sided p-value:           0.0001

Details:
    number of observations:   [10,20]
    t-statistic:              -4.493473545511074
    degrees of freedom:       28
    empirical standard error: 0.29053257978196967

等分散を仮定しないT検定

julia> x1 = rand(Normal(0,3), 10)
julia> x2 = rand(Normal(1,1), 20)
julia> UnequalVarianceTTest(x1, x2)
Two sample t-test (unequal variance)
------------------------------------
Population details:
    parameter of interest:   Mean difference
    value under h_0:         0
    point estimate:          -1.8905345365179094
    95% confidence interval: (-4.01, 0.229)

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.0750

Details:
    number of observations:   [10,20]
    t-statistic:              -1.9826311751666585
    degrees of freedom:       10.182596382128017
    empirical standard error: 0.9535482747359667

3.カイ二乗検定

とりあえずクロス集計表を作る

julia> p1 = 0.01
julia> n1 = 10000
julia> Ber1 = Bernoulli(p1)
julia> x1= rand(Ber1, n1)
julia> A = [sum(x1) n1-sum(x1)]
1×2 Array{Int64,2}:
 100  9900

julia> p2 = 0.015
julia> n2 = 8000
julia> Ber2 = Bernoulli(p2)
julia> x2= rand(Ber2, n2)
julia> B = [sum(x2) n2-sum(x2)]
1×2 Array{Int64,2}:
 138  7862

julia> X = [A ; B]
2×2 Array{Int64,2}:
 100  9900
 138  7862

検定

julia> ChisqTest(X)
Pearson's Chi-square Test
-------------------------
Population details:
    parameter of interest:   Multinomial Probabilities
    value under h_0:         [0.007345679012345679, 0.005876543209876543, 0.5482098765432099, 0.4385679012345679]
    point estimate:          [0.005555555555555556, 0.007666666666666666, 0.55, 0.43677777777777776]
    95% confidence interval: Tuple{Float64,Float64}[(0.0, 0.0132), (0.0, 0.0154), (0.5423, 0.5577), (0.4291, 0.4445)]

Test summary:
    outcome with 95% confidence: reject h_0
    one-sided p-value:           <1e-4

Details:
    Sample size:        18000
    statistic:          17.904808584847267
    degrees of freedom: 1
    residuals:          [-2.802226672320168, 3.132984663835462, 0.3243738295754118, -0.3626609665262763]
    std. residuals:     [-4.231407400008412, 4.231407400008412, 4.2314074000084645, -4.2314074000084645]

他の検定も概ね網羅している模様。
ただし、ANOVA&効果量と検出力も見つからなかった...。

chapter.4 可視化

可視化のライブラリは色々ある模様。
ここでは StatsPlots.jl を使った可視化を行います。
またこの chapter は jupyter notebook 前提で書いています。
詳細:https://docs.juliaplots.org/latest/tutorial/

1.パッケージ install

import Pkg; Pkg.add("StatsPlots")
using StatsPlots
using Distributions

2.ヒストグラム

norm_dist = Normal() 
norm = rand(norm_dist,1000)
StatsPlots.histogram(norm, bins=20)

2つのヒストグラムを重ねる

norm_dist1 = Normal(0)
norm_dist2 = Normal(3)
norm1 = rand(norm_dist1,1000)
norm2 = rand(norm_dist2,1000)
StatsPlots.histogram([norm1 norm2])

p1.png

! をつけると後から追加できる。

norm_dist1 = Normal(0)
norm_dist3 = Normal(5)
norm3 = rand(norm_dist3,1000)
StatsPlots.histogram!(norm3)

p2.png

3.折れ線グラフ

coscurve = cos.(-5.0:0.1:5.0)
StatsPlots.plot(coscurve)

p4.png

2つ重ねる

sincurve = sin.(-5.0:0.1:5.0)
StatsPlots.plot([coscurve sincurve])

p5.png

4.散布図

N = 100
norm_noize = rand(Normal(0,0.1) , N)
x = rand(Normal(0,2) , N)
y = sin.(x) .+ norm_noize
StatsPlots.plot(x,y,seriestype=:scatter)

p6.png

こちらも2つ重ねる

y2 = cos.(x) .+ norm_noize
StatsPlots.plot!(x,y2,seriestype=:scatter)

p7.png

5.Boxプロット

y = rand(100,4) # Four series of 100 points each
StatsPlots.boxplot(["Series 1" "Series 2" "Series 3" "Series 4"],y,leg=false)

p9.png

6.Violinプロット

(私は昆布グラフと呼んでいる。あんまり通じない)

y = rand(100,4)
StatsPlots.violin(["Series 1" "Series 2" "Series 3" "Series 4"],y)

p10.png

7.複数同時に表示

例1.

x = 1:10; y = rand(10,4) 
StatPlots.plot(x, y, layout=(4, 1))

p11.png

例2.

x = 1:10; y = rand(10, 4) 
p1 = StatPlots.plot(x, y)
p2 = StatPlots.scatter(x, y)
p3 = StatPlots.plot(x, y)
p4 = StatPlots.histogram(x, y)
StatPlots.plot(p1, p2, p3, p4, layout=(2, 2), legend=false)

p12.png

※棒グラフわからんかった...。
※散布図行列などの一斉に描画するタイプのものも見つからず...。


以上
次回は データフレーム操作の予定です。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした