2019年5月26日日曜日

[Statistics] 標準偏差

Julia, Python, Rで標準偏差の計算を行う際のメモです.

標準偏差に関しては,下式
$\sigma = \sqrt{\sigma^2} = \sqrt{\cfrac{1}{N} ( x_i - \mu )^2}$
の標本分散の平方根とする場合と,下式
$\sigma^2 = \cfrac{1}{N-1} ( x_i - \mu )^2$
の不偏分散の平方根と説明してある場合があります.
名称に関してはどう呼ぶのが正しいのかはわかりませんが,以下の説明では最初の式を標本の標準偏差,2番目の式を不偏標準偏差と呼ぶことにします.

Julia, Python, Rの各言語の標準偏差を求める関数を用いたときに標本の標準偏差と不偏標準偏差のどちらが返ってくるのかを調べた結果です.結論としては,3種類の言語ともに(デフォルトでは)不偏標準偏差の値を返してきます

Julia(JuliaのStatisticsパッケージのドキュメントはこちら
Statistics パッケージの std 関数で標準偏差を求めると,不偏標準偏差が計算されます.これを標本の標準偏差だと思っていると当然,計算結果が合わないので悩むことになります.

以下に,その例を示します.
まずは,ターミナルからJuliaを立ち上げて以下のように1〜10の平均を求めます.
julia> using Statistics
julia> x = collect(1:1:10)         #[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
10-element Array{Int64,1}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
julia> μ = Statistics.mean(x)     # \mu[tab]
5.5
julia> n = length(x)
10

ここで,定義に従って分散を求めると以下の値になります.
julia> sum((x .- μ).^2) / n
8.25
不偏分散を求めると以下の値となります.
julia> sum((x .- μ).^2) / (n - 1)
9.166666666666666

標本の標準偏差を求めると以下の値になります.
julia> sqrt(sum((x .- μ).^2) / n)

2.8722813232690143


不偏標準偏差を求めると以下のようになります.
julia> sqrt(sum((x .- μ).^2) / (n - 1))

3.0276503540974917


続いて,JuliaのStatisticsモジュールのvar(分散)関数を使うと以下のようになります.
julia> s²₂ = Statistics.var(x)

9.166666666666666

この値は不偏分散です(この例では,Statistics.var としていますが,var()だけでも良い).

std(標準偏差)関数を使うと以下のようになります.
julia> s₂ = Statistics.std(x)
3.0276503540974917
この値は不偏標準偏差の値です(この例では,Statistics.std としていますが,std()だけでも良い).

以上から,JuliaのStatisticsモジュールのstd関数を用いると不偏標準偏差の値が返ってくることが確認できました.Jupyter Notebookでの実行例はGiHubで見ることができます.コードもGitHubで見ることができます.

Python(Numpyのstd関数に関するドキュメントはこちら
Juliaと同様にPythonでも確認してみます.
>>> import numpy as np

標準偏差を求める配列を設定します.
>>> x = [i for i in range(1,11)]

デフォルトでのstd関数を用いた場合の値を求めます.
>>> s2 = np.std(x) # np.std(x, ddof=0)
>>> print ("s2 = ", s2)
s2 =  2.8722813232690143
上記のs2の計算式をしている行でもコメントアウトしていますが,デフォルトではddof = 0として計算を行っており,不偏標準偏差の値が返ってきています.

標本の標準偏差を求める場合には,以下のように引数の設定でddof = 1とします.
>>> s1 = np.std(x, ddof=1)
>>> print ("s1 = ", s1)
s1 =  3.0276503540974917

以上から,PythonのNumpyモジュールのstd関数を用いると不偏標準偏差の値が帰ってくることが確認できました.Jupyter Notebookでの実行例はGitHubで見ることができます.コードもGitHubで見ることができます.

R
Rでも確認してみます.
まずは,標準偏差を求める配列を設定します.
> x = (1:10)

標準偏差を求めてみます.
> s2 = sd(x)
> s2
[1] 3.02765
不偏標準偏差の値が返ってきています.

標本の標準偏差を求めるには以下のようにします.
> v = var(x)*(length(x)-1)/length(x)
> s1 = sqrt(v)
> s1
[1] 2.872281

以上から,Rのsd関数を用いると不偏標準偏差の値が帰ってくることが確認できました.Jupyter Notebookでの実行例はGitHubで見ることができます.コードもGitHubで見ることができます.

0 件のコメント :

コメントを投稿