Colorless Green Ideas

R言語で母集団の分散・標準偏差を求める方法

概要
R言語に標準で付属しているパッケージでは、母集団と見なしたときの分散・標準偏差を求める関数がない。一方で、GE / sjstats / smallstuff といったパッケージを使えば、母集団の分散・標準偏差を求めることができる。

はじめに

R言語に標準で付属しているパッケージの1つ stats には分散・標準偏差を求める関数として var() と sd() がある。これらは標本における不偏分散・不偏標準偏差を返す関数である。要するに、$n – 1$ で割る形の分散・標準偏差を求める関数である。

一方で、標準で付属しているパッケージには、母集団における分散・標準偏差（$n$ で割る形のもの）を求める関数はない。このため、母集団の分散・標準偏差を計算できる関数を含むパッケージを使うか、自作関数を用意しておく必要がある。

以下では、母集団における分散・標準偏差を計算できるパッケージを3つ紹介する（GE / sjstats / smallstuff）。

GE パッケージ

GE パッケージには、母集団の分散・標準偏差を求める関数として、var.p() と sd.p() が用意されている。

このパッケージは CRAN に掲載されており、以下のコマンドでインストールできる。

install.packages("GE")

以下の例では、GE パッケージを読みこんで、数値が入ったベクトル x を母集団とみて分散・標準偏差を計算している。

library("GE")

x <- c(3, 5, 7, 12)
var.p(x) # [1] 11.1875
sd.p(x)  # [1] 3.344772

なお、欠損値を含む数値を与えた場合、これらの関数は、欠損値 NA を返す。欠損値を除外し、残った値だけで計算したければ、オプションで na.rm = TRUE と指定すればよい。))。これは標本における不偏分散・不偏標準偏差を返す var() と sd() と同じ挙動である。

x <- c(1, 4, 12, NA, 23)

var.p(x) # [1] NA
sd.p(x)  # [1] NA

var.p(x, na.rm = TRUE) # [1] 72.5
sd.p(x, na.rm = TRUE)  # [1] 8.514693

さらに、var.p() と pop.sd.p() は、wt オプションで指定することで、重み付け有りの計算ができる。以下の例では3番目の数値 (7) の重みを他の数値の2倍にして計算している。

x <- c(3, 5, 7, 12)
w <- c(1, 1, 2, 1)

var.p(x, wt = w) # [1] 8.96
sd.p(x, wt = w)  # [1] 2.993326

smallstuff パッケージ

smallstuff パッケージには、母集団の分散・標準偏差を求める関数として、pop.var() と pop.sd() が用意されている。なお、これらの関数は、欠損値を自動的に除外して計算する。

このパッケージも CRAN に掲載されており、以下のコマンドでインストールできる。

install.packages("smallstuff")

以下の例では、smallstuff パッケージを読みこんで、数値が入ったベクトル x を母集団とみて分散・標準偏差を計算している。

ibrary("smallstuff")

x <- c(1, 4, 12, NA, 23)
pop.var(x) # [1] 72.5
pop.sd(x)  # [1] 8.514693

sjstats パッケージ

sjstats には、母集団の分散・標準偏差を求める関数として、var_pop(x) と sd_pop(x) が用意されている。

このパッケージも CRAN に掲載されており、以下のコマンドでインストールできる。

install.packages("sjstats")

以下の例では、sjstats パッケージを読みこんで、数値が入ったベクトル x を母集団とみて分散・標準偏差を計算している。マニュアルには明記されていないが、これらの関数も欠損値を自動的に除外して計算する。

library("sjstats")

x <- c(1, 4, 12, NA, 23)
var_pop(x) # [1] 72.5
sd_pop(x)  # [1] 8.514693

各パッケージの関数の定義の違い

上で紹介した3つのパッケージでの分散の関数の定義はそれぞれ少しずつ異なっている。

GE パッケージでは、 $\frac{\sum_{i=1}^{n} (x_i – \mu)^2}{n}$ となるように計算している。つまり、各データが母平均からどれくらい離れているか（偏差）を二乗したものの平均を取っている。これは最も一般的な分散の定義に従って計算している。Rの式で書くと、mean(sum(x^2) - mean(x)^2) のようになる。

smallstuff パッケージでは、$\frac{\sum_{i=1}^{n} x_i^2 – n \mu^2}{n}$ となるように定義している。要するに、平方和（各値の二乗の総和）から、母平均の平方に母集団の大きさ $n$ をかけたものを引き、その結果を $n$ で割っている。Rの式で書くと、(sum(x^2) - length(x) * mean(x)^2)/ length(x) のようになる。

sjstats パッケージでは、単に stats::var(x) * ( (n - 1) / n) のように計算している。つまり、標本における不偏分散に $\frac{n-1}{n}$ をかけている。要するに、元々 $n-1$ で割ったものを打ち消し、$n$ で割り直しているということだ。ただ、var() を使っている副作用として、要素が1つしかない場合に、分散・標準偏差として NA が返ってきてしまう。

これらは、数学的に言えば同じ値を返すはずだ。だが、実際の計算の際には丸めのタイミングなどが違うために、異なった値となってしまうことがある。例えば、（あくまで私の環境での例だが）以下では小数第14位以下の数字が異なる結果となった。

x <- c(0.000001, 4.0002, 12, 25)
sprintf("%.40g", GE::var.p(x))           # [1] "91.186869882475178883"
sprintf("%.40g", sjstats::var_pop(x))    # [1] "91.186869882475193094"
sprintf("%.40g", smallstuff::pop.var(x)) # [1] "91.186869882475164673"

このように計算結果が異なることは、自分の計算結果と他の人の計算結果を突き合わせる際に問題になりうる。例えば、計算ミスがないように自分と他のもう1人で同じデータに対して別々に計算を進め、最後に両者が合っているか確認するものとする。この場合、双方が別々のパッケージを用いていると、同じデータでも結果が異なってしまうということがありうるのだ。

逆に言うと、結果がわずかに異なってしまう場合、数学的には同等な別の定義で計算している可能性を疑って調べてみると解決することも少なくない。これはRの異なるパッケージ同士の話だが、他の統計ソフトとRを比べるときにもこういうことはありうる。

注釈