Julia で統計学分布関数と仮説検定

(@v1.5) pkg> add Distributions, Random, HypothesisTests
julia> using Distributions, Random, HypothesisTests

分布関数

標準正規分布において，R の rnorm, pnorm, qnorm, dnorm に相当するもの

Julia の場合は，「オブジェクト」を操作する。

まず，標準正規分布のオブジェクトを作る。

julia> obj = Normal()
Normal{Float64}(μ=0.0, σ=1.0)

乱数生成

正規乱数 500 個を生成する。

R では

> d = rnorm(500)

Julia では

julia> d = rand(obj, 500);

生成された乱数データについて，平均値，不偏分散，標準偏差を計算してみる。

julia> println("mean = $(mean(d)), var = $(var(d)), sd = $(std(d))")
mean = -0.025412744391597998, var = 0.9232796228593824, sd = 0.9608744053513875

下側確率

R では

> pnorm(1.96)
[1] 0.9750021

Julia では

julia> cdf(obj, 1.96)
0.9750021048517795

複数の引数（ベクトル）において，各要素についての結果は，ドット記法 'cdf.(' を使って書ける。

R では

> pnorm(c(-1.96, 0, 1.96))
[1] 0.0249979 0.5000000 0.9750021

Julia では

julia> cdf.(obj, [-1.96, 0, 1.96])
3-element Array{Float64,1}:
 0.024997895148220435
 0.5
 0.9750021048517795

パーセント点

R では

> qnorm(c(0.025, 0.5, 0.975))
[1] -1.959964  0.000000  1.959964

Julia では

julia> quantile(obj, [0.025, 0.5, 0.975])
3-element Array{Float64,1}:
 -1.9599639845400592
  0.0
  1.9599639845400583

確率密度

R では

> dnorm(c(-1.96, 0, 1.96))
[1] 0.05844094 0.39894228 0.05844094

Julia では

julia> pdf.(obj, [-1.96, 0, 1.96])
3-element Array{Float64,1}:
 0.05844094433345147
 0.3989422804014327
 0.05844094433345147

等分散を仮定した独立 2 標本の平均値の差の検定

両側検定

R では

> x = c(3,2,1,2,3,2,1,2,3,4,3,2,2,3,4)
> y = c(3,2,1,2,3,2,3,4,4,3,3,3,2,2,3,3,4,5)
> t.test(x, y,  var.equal=TRUE)

 Two Sample t-test

data:  x and y
t = -1.282, df = 31, p-value = 0.2093
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.0939215  0.2494771
sample estimates:
mean of x mean of y 
 2.466667  2.888889 

Julia では

julia> x = [3,2,1,2,3,2,1,2,3,4,3,2,2,3,4];
julia> y = [3,2,1,2,3,2,3,4,4,3,3,3,2,2,3,3,4,5];
julia> obj = EqualVarianceTTest(x, y)
Two sample t-test (equal variance)
----------------------------------
Population details:
    parameter of interest:   Mean difference
    value under h_0:         0
    point estimate:          -0.422222
    95% confidence interval: (-1.0939, 0.2495)

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

Details:
    number of observations:   [15,18]
    t-statistic:              -1.2820140381935845
    degrees of freedom:       31
    empirical standard error: 0.32934290081343587

片側検定

R で，alternative="less"

> t.test(x, y,  var.equal=TRUE, alternative="less")

 Two Sample t-test

data:  x and y
t = -1.282, df = 31, p-value = 0.1047
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
      -Inf 0.1361849
sample estimates:
mean of x mean of y 
 2.466667  2.888889 

Julia では

julia> println(pvalue(obj, tail=:left))
0.10467146197969886

R で，alternative="greater"

> t.test(x, y,  var.equal=TRUE, alternative="greater")

 Two Sample t-test

data:  x and y
t = -1.282, df = 31, p-value = 0.8953
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
 -0.9806293        Inf
sample estimates:
mean of x mean of y 
 2.466667  2.888889 

Julia では

julia> println(pvalue(obj, tail=:right))
0.8953285380203011

その他の検定手法

OneSampleTTest(x)
EqualVarianceTTest(x, y)
UnequalVarianceTTest(x, y)
ChisqTest(matrix)

これ以外は，GitHub にあったりするが，最悪，自分で書かねばならない。かな？