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

最近流行りつつある科学計算処理向けの高水準言語「Julia」の基本操作をまとめていきます。
1.[準備編](https://qiita.com/HiroyukiTachikawa/items/f4c7a620987a45b0429b)
2.[確率分布・仮説検定と可視化編](https://qiita.com/HiroyukiTachikawa/items/dfc1b0e92d80287ea7f6)
3.[データフレーム操作編](https://qiita.com/HiroyukiTachikawa/items/e01917ade931031ec6a1)
## chapter.2 確率分布
確率分布のパッケージがあったので、ほぼそのパッケージの紹介のような感じです。
R とはまた使用感が異なりますが、個人的にはこっちの方が直感的で好きです。

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

### 1.パッケージ install
```julia
julia> import Pkg; Pkg.add("Distributions")
julia> using Distributions
```
この後、ランダムにデータを生成する場面がありますが、seed を固定したい場合は以下

```julia
julia> using Random
julia> Random.seed!([シード値])
```
### 2.確率分布オブジェクト
とりあえず名称そのま関数を実行してみるとわかる通り、デフォルトのパラメーターが設定されている。もちろん指定することもできる。

```julia
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
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
julia> println("μ : ", norm_dist.μ)
μ : 50.0

julia> println("σ : ", norm_dist.σ)
σ : 10.0
```
配列の場合

```julia
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
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
julia> norm_fitted = fit(Normal, array) # ここでfit
Normal{Float64}(μ=-0.5830417798446266, σ=0.8450652913045483)
```

確認してみる

```julia
julia> println(norm_fitted)
Normal{Float64}(μ=-0.5830417798446266, σ=0.8450652913045483)
```
パラメータを取り出す場合は先ほど同様...

```julia
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
julia> import Pkg; Pkg.add("HypothesisTests")
julia> using HypothesisTests, Distributions
```
### 2.とりあえずT検定
一標本

```julia
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
julia> println("right-side p-value : ", pvalue(ttest_object, tail=:right))
right-side p-value : 0.0005084438096162109
```
一応 Paired T

```julia
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
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
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
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
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]
```
他の検定も概ね網羅している模様。

## chapter.4 可視化
可視化のライブラリは色々ある模様。
ここでは StatsPlots.jl を使った可視化を行います。
またこの chapter は jupyter notebook 前提で書いています。
詳細:https://docs.juliaplots.org/latest/tutorial/
### 1.パッケージ install
```julia
import Pkg; Pkg.add("StatsPlots")
using StatsPlots
using Distributions
```
### 2.ヒストグラム

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

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

```julia
norm_dist1 = Normal(0)
norm_dist2 = Normal(3)
norm1 = rand(norm_dist1,1000)
norm2 = rand(norm_dist2,1000)
StatsPlots.histogram([norm1 norm2])
```
![p1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/112604/737874d4-8230-b7ae-e7f7-07235feec776.png)

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

```julia
norm_dist1 = Normal(0)
norm_dist3 = Normal(5)
norm3 = rand(norm_dist3,1000)
StatsPlots.histogram!(norm3)
```
![p2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/112604/d01ac4ff-e117-400f-ecb8-d76bab270a11.png)

### 3.折れ線グラフ
```julia
coscurve = cos.(-5.0:0.1:5.0)
StatsPlots.plot(coscurve)
```
![p4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/112604/51fbc712-d964-600e-b4ff-311ae9ba5bd5.png)

2つ重ねる

```julia
sincurve = sin.(-5.0:0.1:5.0)
StatsPlots.plot([coscurve sincurve])
```
![p5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/112604/baa04d53-ed49-5882-fc6d-2521cbb9d49e.png)

### 4.散布図
```julia
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/112604/f144fb5b-8761-98c7-bebd-17c3a2f9c1ec.png)

こちらも2つ重ねる

```julia
y2 = cos.(x) .+ norm_noize
StatsPlots.plot!(x,y2,seriestype=:scatter)
```
![p7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/112604/2be3e8cb-080c-3a81-8c41-271f4bd97f98.png)

### 5.Boxプロット
```julia
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/112604/4192e537-87ee-bcd4-4ee4-119ba8313408.png)

### 6.Violinプロット
(私は昆布グラフと呼んでいる。あんまり通じない）

```julia
y = rand(100,4)
StatsPlots.violin(["Series 1" "Series 2" "Series 3" "Series 4"],y)
```
![p10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/112604/d2a892d8-110f-41d9-39a7-c57ad640bd19.png)

### 7.複数同時に表示

例1.

```julia
x = 1:10; y = rand(10,4) 
StatPlots.plot(x, y, layout=(4, 1))
```
![p11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/112604/e54dde93-049c-b92d-4a0b-956f0386c4af.png)

例2.

```julia
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/112604/cfffd3ca-ec1e-476b-6d58-98eba063edaf.png)

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

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