目次

・私たちは母集団を知りたい

・基本的な統計量のおさらい

　　・標本

　　・標本平均、標本分散

　　・不偏分散、標準偏差

・不偏推定量の妥当性をシミュレーションする

・なぜ標本分散は過少推定となるか？

・不偏分散の平方根（標準偏差）は不偏推定量ではない！！

・自分が何を計算しているのか把握しよう！

・補足1：標本分散が不偏推定量ではないことを証明

・補足2：外れ値の影響を受けにくい統計量

　　・中央値、四分位数

　　・最頻値

私たちは母集団を知りたい

　私たちが解析をするときに、平均や分散などの統計量を求めたり、検定を行うのはなぜだろうか？　それはずばり「標本の母集団・真の状態を知りたい」と考えているからだ。手元の標本は真の状態の一部にしか過ぎない。しかし、真の状態を知るために、母集団すべてを調べることは不可能だ。そこで、私たちは標本から得られる統計量を調べることで、真の状態を推定しようとしている。この過程が、平均や分散などの統計量を求めたり、検定を行うことなのだ。

直感的には、標本は母集団の一部なのだから、母集団の代表値が取れていれば標本から計算した平均（標本平均）や分散（標本分散）は、母集団の平均（母平均）や分散（母分散）と一致するだろう。というか、一致していてほしい。ゆえに、標本平均や標本分散は母集団を知るための手がかりとして活用される。しかし、統計を学んでいると、分散に関しては不偏分散なるものも習う。計算方法もよく似ていて非常にややこしい。そこで本ページでは、統計量は母集団の状態を推定する道具だという点に着目して、平均や分散の性質を紹介する。結論から言えば、標本分散よりも不偏分散のほうが、母分散のよりよい推定値になっているのだ。

基本的な統計量のおさらい

・標本（サンプル）

まず、標本（サンプル）という言葉を押さえよう。統計学における1つの標本とは、同じ母集団に由来すると思われる観測値の集まりだ。検定などで比較対象となる場合は、標本のことを群、あるいは分散分析や線形モデルでは要因と呼ぶこともある。勘違いされやすいが、1つの標本=1つの観測値ではない。ゆえに標本数というと「観測値の集まりの個数」を指すのであって、標本の中の観測値の個数を指すわけではない。

　「標本の中の観測値の個数」は標本の大きさ（サンプルサイズ）と呼ばれる。上記のような使われ方をしない場面もあるが、統計学で考えるときは上記のように定義しよう。そうしないと、以降の議論がすべて破綻してしまう。

・標本平均、標本分散

標本の要約に用いられる代表的な統計量は標本平均と標本分散だ（ほかの統計量は追記する）。ともに標本の母集団の母平均と母分散の推定値となる。標本平均は分布の「中心」の指標値であり、標本分散は分布の幅広さや標本のばらつきの指標値とされる。

なぜ分散は上記のように定義されるのだろう。いいかえれば、分布の幅広さとか標本のばらつきの指標として、適切だと言えるだろうか。それは、以下の図を見るとよいだろう。

例え、平均が等しいとしても、分散が違えば、標本の散らばり方=分布の幅広さが異なることが一目瞭然だ。分散が大きければ大きいほど、分布は幅広くなったり、極端な値が多くなったりすることがわかる。

・不偏分散、標準偏差

不偏分散も、母分散の推定値だ。また、不偏分散の平方根をとったものを標準偏差と呼ぶ。

標本分散と不偏分散を見比べると、偏差の平方和を割る数が異なっている。標本分散はサンプルサイズのままだが、不偏分散ではサンプルサイズから1引いた値で割っているのだ。この微妙な違いと散見される用語の誤用により、統計初学者は混乱する。平均はnで割っているし、分散もnで割りたくなる気持ちはわかる。用語も、n-1で割る不偏分散のことを標本分散と呼んでいることがある。こんなことがあると、ますます標本分散と不偏分散は何が違うのか、わからなくなってしまう。ともかく、このページではnで割っているほうを標本分散、n-1で割っているほうを不偏分散と呼ぶ。

　では、この2つは何が違うのかというと、標本の分散を計算するにあたっては標本分散よりも不偏分散のほうが母分散に近い値を示すのだ。つまり、不偏分散のほうが母分散の推定値としては適切であるといえる。ただし、母分散を推定する必要がない場合、例えば、「標本のばらつきを知りたいだけ」とか、「標本と母集団が一致しているとき」は、標本分散を使って差し支えはない。母分散を推定したいときに、不偏分散を使うのだ。解析の記録を残すときは、自分がどういう目的で、どういう計算をしたのか、を明確にしておく必要があるので、心しておこう。

不偏推定量の妥当性をシミュレーションする

　まず、標本平均が母平均の推定値として妥当か、をRをつかったシミュレーションで明らかにしていく。そのあとに、標本分散よりも不偏分散のほうが母分散に近い値を示す、ということを確かめよう。

　標本平均や不偏分散のように、サンプルサイズに関わらず、その期待値が母集団の統計量と一致する統計量のことを不偏推定量 unbiased estimatorと呼ぶ。不偏分散の「不偏」はここからきている。標本平均も、普通はそう呼ばないが、不偏平均と呼ぶこともできるだろう。一方で、標本分散は、母集団の統計量とは一致せず、やや母分散よりも小さな値となる。つまり標本分散は、小さくなる偏りbiasがあり、その点で不偏ではないのだ。しかし、サンプルサイズが大きくなると母集団の統計量に収束する。このような統計量は一致推定量 consistent estimatorと呼ぶ。収束するというのはつまり、サンプルサイズが大きければ、標本分散も不偏分散も値が大きく変わらなくなるということだ。例えば、サンプルサイズn = 2だと、1で割るか、2で割るかは大きな違いとなるが、n = 10000なら、9999で割ろうが、10000で割ろうが大差はないだろう。

　前置きが長くなったが、標本平均のシミュレーションを行おう。まずは、母集団を定義する。今回は、平均が0、標準偏差が1の正規分布だ。

------------------------------------------------------

library(plotn)

#このライブラリは私オリジナルの作図ライブラリである。

#もし、plotnライブラリが欲しい場合は、

#以下のようにコマンドを打つ（#は消す）。devtoolsライブラリが必要。

#library(devtools)

#install_github("bugplant/plotn")

n <- rnorm(10000, mean = 0, sd = 1) #平均0、標準偏差1の正規分布に従う乱数を10000個生成

histn(n, xlab = "value", ylab = "頻度") #作図。シミュレーションなので乱数を生成するごとに少しずつ図は異なってくる。hist(n)でも可。

------------------------------------------------------

このような母集団からサンプルサイズn = 3の標本が生成され、その標本平均を算出する過程を10000回繰り返すとしよう。そのコードは以下の通りだ。

------------------------------------------------------

m <- NULL

for(i in 1:10000) {

  x <- rnorm(3, mean = 0, sd = 1) #平均0、標準偏差1の正規分布からサンプルサイズ3の標本を生成

  m <- c(m, mean(x))#平均の計算

}

histn(m, xlab = "m", ylab = "頻度", xaxt = "n")#標本平均の分布を図示

m_mean <- mean(m)#標本平均の平均の計算

overdraw(abline(v = m_mean, col = "red"), #標本平均の平均を図示

         axis(side = 1, at = -5:5, labels = -5:5, cex.axis = 1.1))

m_mean #標本平均の平均の表示。ほぼ0となっていることを確認

## [1] 0.003011211

------------------------------------------------------

確かに図2における赤線、すなわち標本平均を10000回算出した時のその平均はほとんど0となっている。小さなサンプルサイズであることの影響を受けているようには見えない。これを見れば標本平均が、「平均的に見れば」母平均に近くなる、つまり不偏推定量であるということを確認できるのではないだろうか。

　続いて、同じ手続きを標本分散についても行ってみよう。母集団の標準偏差は1なので、その2乗の1が母分散となる。もし標本分散が不偏推定量であるなら、上記と同様に、10000回標本分散を算出した時の、標本分散の平均は1となるはずである。

------------------------------------------------------

s_2 <- NULL

for(i in 1:10000) {

  x <- rnorm(3, mean = 0, sd = 1)#平均0、標準偏差1の正規分布からサンプルサイズ3の標本を生成

  s_2 <- c(s_2, 1/3 * sum((x - mean(x))^2))#標本分散の計算。割る数がサンプルサイズと同じであることに注目

}

histn(s_2, xlab = "s^2", ylab = "頻度", xaxt = "n")#標本分散の分布を図示

s2_mean <- mean(s_2)#標本分散の平均の計算

overdraw(abline(v = s2_mean, col = "red"),#標本分散の平均を図示

         axis(side = 1, at = 0:10, labels = 0:10, cex.axis = 1.1))

s2_mean #標本分散の平均の表示。母分散の2/3倍となっていることを確認

## [1] 0.6614964

------------------------------------------------------

ところが標本分散の平均（赤線）は、母分散とは一致せず、より小さく推定される。今回の例だと、大体0.66であり、母分散の約2/3倍程度だ。この2/3という値は、標本分散と不偏分散の比率、(n-1)/nに由来している。では、不偏分散を同様に計算して、不偏推定量であることを確認しよう。

------------------------------------------------------

u_2 <- NULL

for(i in 1:10000) {

  x <- rnorm(3, mean = 0, sd = 1) #平均0、標準偏差1の正規分布からサンプルサイズ3の標本を生成

  u_2 <- c(u_2, 1/2 * sum((x - mean(x))^2)) #不偏分散の計算。割る数がサンプルサイズから1引いた数であることに注目

}

histn(u_2, xlab = "u^2", ylab = "頻度", xaxt = "n") #不偏分散の分布を図示

u2_mean <- mean(u_2)　#不偏分散の平均の計算

overdraw(abline(v = u2_mean, col = "red"), #不偏分散の平均を図示

         axis(side = 1, at = 0:10, labels = 0:10, cex.axis = 1.1))

u2_mean #不偏分散の平均の表示。ほぼ1となっていることを確認

## [1] 0.9941563

------------------------------------------------------

確かに、不偏分散の平均（赤線）は母分散と一致する。このように標本分散と不偏分散はよく似た統計量であるが、推定の性質が異なることは押さえておこう。

なぜ標本分散は過少推定となるか？

　次に標本分散が過少推定になってしまう原因を考察することで、推定に関する理解を深めていこう。一言で言えばその原因は、標本平均が母平均からずれてしまうためだ（数値的証明は補足）。図2を見るとわかるが、標本平均の平均は母平均に一致するし、標本平均は母平均を中心とした釣り鐘型分布となって母平均に近い値を示すことが多い（中心極限定理）。しかし、母平均±1くらいの範囲なら、そこそこの確率でずれる可能性がある。このずれが、標本分散の過少推定を引き起こす。これは、標本分散の定義を微分することで確認できる。

標本分散の定義式は一見複雑だが、微分して確認してみると、標本平均を頂点とする下に凸の2次関数だとわかる。つまり、常に（標本平均を使った標本分散）≦（母平均を使った"標本分散"）となるのだ。私たちは母平均を直接知ることができないから、標本平均を使って推定しようとしている。したがって、標本平均を使うことで標本分散が母分散からずれてしまうことは避けられない。だからこそ、正しい推定のために不偏分散のような理論から導かれた補正値が必要になるわけだ。

　もし、私たちが母平均を知っていたらどうなるだろう。上記の議論から考えると、母平均を使って”標本分散”を計算すれば、母分散と一致すると予測できる。まさに、今のシミュレーションの状況がそんな状況だ。私たちは母集団を定義しているわけだから母平均も母分散も分かった状態で分析をしている。

------------------------------------------------------

s_2 <- NULL

for(i in 1:10000) {

  x <- rnorm(3, mean = 0, sd = 1)#平均0、標準偏差1の正規分布からサンプルサイズ3の標本を生成

  s_2 <- c(s_2, 1/3 * sum((x - 0)^2))#母平均0を使って"標本分散"を計算。

}

histn(s_2, xlab = "s^2", ylab = "頻度", xaxt = "n")#"標本分散"の分布を図示

s2_mean <- mean(s_2)#"標本分散"の平均の計算

overdraw(abline(v = s2_mean, col = "red"),#"標本分散"の平均を図示

         axis(side = 1, at = 0:10, labels = 0:10, cex.axis = 1.1))

s2_mean #"標本分散"の平均の表示。ほぼ1となっていることを確認

## [1] 0.993571

------------------------------------------------------

予想通り、母平均を使った"標本分散"の平均（赤線）は母分散と一致する。しかし、繰り返しになるが、実用上は母平均を知ることはできない。ゆえに、不偏分散を使って母分散を推定することになる。

不偏分散の平方根（標準偏差）は不偏推定量ではない！！

　ここからは、ちょっと発展的な内容となる。あまり議論の対象となることがないが、不偏推定量を学ぶ上で示唆に富んでいると感じた話題だ。このページの最初のほうで標準偏差を不偏分散の平方根として定義した。ほかに2通り、標準偏差の定義がある。その1つが、標本分散の平方根を標準偏差とするものだ。では、不偏分散の平方根と標本分散の平方根のどちらが、母標準偏差の不偏推定量なのだろうか？　実はどちらも不偏推定量ではない。一見すれば、不偏分散が不偏推定量なので、その平方根も不偏推定量となる気がする。確認してみよう。

------------------------------------------------------

s <- NULL

for(i in 1:10000) {

  x <- rnorm(3, mean = 0, sd = 1)#平均0、標準偏差1の正規分布からサンプルサイズ3の標本を生成

  s <- c(s, sqrt(1/3 * sum((x - mean(x))^2)))#標本分散の平方根を計算。

}

histn(s, xlab = "s", ylab = "頻度", xaxt = "n")#標本分散の平方根の分布を図示

s_mean <- mean(s)#標本分散の平方根の平均の計算

overdraw(abline(v = s_mean, col = "red"),#標本分散の平方根の平均を図示

         axis(side = 1, at = 0:10, labels = 0:10, cex.axis = 1.1))

s_mean #標本分散の平方根の平均の表示。母標準偏差の過少推定となっている。

## [1] 0.7239629

u <- NULL

for(i in 1:10000) {

  x <- rnorm(3, mean = 0, sd = 1)#平均0、標準偏差1の正規分布からサンプルサイズ3の標本を生成

  u <- c(u, sqrt(1/2 * sum((x - mean(x))^2)))#不偏分散の平方根を計算。

}

histn(s, xlab = "u", ylab = "頻度", xaxt = "n")#不偏分散の平方根の分布を図示

u_mean <- mean(u)#不偏分散の平方根の平均の計算

overdraw(abline(v = u_mean, col = "red"),#不偏分散の平方根の平均を図示

         axis(side = 1, at = 0:10, labels = 0:10, cex.axis = 1.1))

u_mean #不偏分散の平方根の平均の表示。母標準偏差の過少推定となっている。

## [1] 0.8870197

------------------------------------------------------

なんと不偏分散の平方根は母標準偏差1の過少推定となる。標本分散の平方根よりかは、不偏分散の平方根のほうが推定精度が良いが、不偏分散の平方根は母標準偏差の不偏推定量ではないのだ。ここで3つめの標準偏差の定義を紹介しよう。母標準偏差の不偏推定量、不偏標準偏差と呼ばれるべきこの値は、以下のように計算する（これは、不偏分散が従うカイ二乗分布から算出される結果である）。

上記の式のように不偏標準偏差は不偏分散の平方根に補正の係数がかかった形となっている。この補正の係数が今までとは毛色が違う複雑な形だ。そもそもガンマ関数ってなんだ？と思う人が多いだろう。詳しくは述べないが、ガンマ関数は階乗の計算の一般化ともいわれるものだ。階乗といえば、普通、3! = 3×2×1のように計算し、自然数で定義されるものだが、これを実数で計算できるように改良したものだ。つまり、(1/2)!みたいなものも計算できる。ただし、積分の式で定義されているため、解説には相応の努力が必要であり、省略させてほしい。（こちらで解説）とりあえず、シミュレーションする分には、関数の名前さえ分かっていればよいので、以下にコードを示す。

------------------------------------------------------

D <- NULL

for(i in 1:10000) {

  x <- rnorm(3, mean = 0, sd = 1)

  D <- c(D, sqrt((3-1)/2)*gamma((3-1)/2)/gamma(3/2)*sqrt(1/2 * sum((x - mean(x))^2)))

}

#上記のgamma()部分がガンマ関数である。

histn(D, xlab = "D", ylab = "頻度", xaxt = "n") #不偏標準偏差の分布を図示

D_mean <- mean(D) #不偏標準偏差の平均の計算

overdraw(abline(v = D_mean, col = "red"), #不偏標準偏差の平均を図示

         axis(side = 1, at = 0:10, labels = 0:10, cex.axis = 1.1))

D_mean  #不偏標準偏差の平均の表示。ほぼ1であることを確認

## [1] 0.99998

------------------------------------------------------

確かに不偏標準偏差は母標準偏差と一致する。ゆえに不偏推定量と言って問題なさそうだ。標準偏差は、定義が多い分、分散以上に混乱に満ち溢れているのだ。

自分が何を計算しているのか把握しよう！

　本ページの最初に述べたことを再掲するが、「解析の記録を残すときは、自分がどういう目的で、どういう計算をしたのか、を明確にしておく必要がある」。手計算するときは、式を書き下すだろうから問題ない。問題があるのはプログラムを使うときだ。例えば、Ｒで標準偏差を計算するsd()というコマンドがある。これが、上記の不偏分散の平方根なのか、標本分散の平方根なのか、不偏標準偏差なのか、ちゃんと把握しているだろうか？　Excelにも標準偏差を計算するSTDEV.PやSTDEV.Sがあるが、それぞれが何を計算しているか把握しているだろうか？　わからないのであれば今一度確認しておこう。そして、最低限、どのような計算をしたか、プログラムは何を使ったかなど、記録を怠らないことだ。

補足1：標本分散が不偏推定量ではないことを証明

　標本分散の期待値が母分散よりも小さくなることを示そう。各データ点が独立であることを仮定すれば、下記のように導くことができる。

式変形途中の式(1)がわかりやすいが、式(1)の左項はその期待値が母分散になる。右項は母分散からいくらかの値を引いているが、標本分散が母分散と一致するためには標本平均と母平均が一致していなければならない。実際、常に一致するとは限らないので、平均的には標本平均は母分散より小さいのである。この発想は、上記で紹介した、標本分散を微分した時と類似している。

補足2：外れ値の影響を受けにくい統計量

　標本平均は、標本の外れ値に強く影響を受ける。例えば、下の例のように、ほとんどの観測値は等しいのに、ある一つの観測値の存在によって、劇的に標本平均が変わってしまう。

------------------------------------------------------

d1 <- c(0,1,1,2,2,2,3,3,3,3,3,3,4,4,4,4,4,5,5,5,5,6,6,7,8)

d2 <- c(0,1,1,2,2,2,3,3,3,3,3,3,4,4,4,4,4,5,5,5,5,6,6,7,100)

mean(d1)

## [1] 3.72

mean(d2)

## [1] 7.4

------------------------------------------------------

　標本平均だけを見れば、あたかもd2は全体のデータがd1よりも2倍程度大きいかのような印象を受けるが、その差はたった一つのデータによってもたらされたものだ。そう考えると、ほとんどのデータの分布が同じであるd1とd2を全く別の分布の”中心”を持つ、というのはやや直観に反する。

　そこで、そのような外れ値に影響を受けにくい分布の中心、あるいは代表値を求める方法が知られている。それが、中央値medianと最頻値modeだ。

・中央値median、四分位数quantile

　以下のようにデータを大きさ順に並べる。そのとき、データの量を50%に分ける値を中央値とする。ともに、中央値を境に、m個のデータが存在するようになっている。

　この中央値の考え方をさらに発展させて、特にデータを25:75、50:50、75:25に分ける点を総称して、四分位数quantileと呼ぶ。データを4つに分けるから「四分」位数である。値が小さい方から第一四分位数、第二四分位数＝中央値、第三四分位数と呼び、第三四分位数 - 第一四分位数を四分位範囲と呼ぶ。

中央値も四分位数も、データの生の値ではなく順位に基づいて決められる値であり、それゆえにある一つの値が極端であってもその影響は受けにくい。

　四分位数は箱ひげ図box plotを作るときに用いられるため、極めて重要である。四分位数を見ればある程度データの分布がわかるのが利点だ。

・最頻値mode

　離散値データの場合、標本中、最も個数の多い値を最頻値とする。

しかしながら、連続値データの場合、厳密には全く同じ値をとるデータは出現しない。そこで、ヒストグラムなどを書いて、もっとも出現頻度の高い領域のデータを最頻値としたり、データからカーネル密度推定することで、もっとも確率密度の高い値を最頻値としたりする。

以上のような、中央値、最頻値、そして、平均値は、左右対称なデータであれば、どれも似たり寄ったりの値になる。しかし、例えば、右に裾が長ければ最頻値<中央値<平均値となり、逆に左に裾が長ければ逆の関係になる。

　以下のようなデータで確認してみよう。ガンマ分布で右に裾の長いデータを生成してみた。

------------------------------------------------------

d <- rgamma(n = 10000, shape = 4, scale = 100)

histn(d, kernel = T)#図6の描画

overdraw(abline(v = mean(d), col = "red"),

         abline(v = median(d), col = "blue"),

         abline(v = k$x[which.max(k$y)], col = "green"))

------------------------------------------------------

確かに、右に裾が長ければ、最頻値<中央値<平均値となっている。具体的には、以下のように、平均値は400、中央値は368、最頻値は302程度だとわかる。

------------------------------------------------------

mean(d)#平均値

## [1] 400.748

median(d)#中央値

## [1] 368.1376

k <- density(d)#カーネル密度推定

k$x[which.max(k$y)]#最頻値

## [1] 301.8384

------------------------------------------------------

よりリアリティーのあるデータとしてとらえるなら、例えばある会社に勤めたときの年収とかどうだろう（0付近がいるのやばすぎだけど）。ある会社が売り文句として、うちの平均年収は400万ですよー、と就活生にアピールする。まあ、400万なら悪くないかと、入社を決める。しかし、年収の分布は図6のようならどうだろう。つまり、会社員の半分は400万に届くどころか、368万円以下の年収で、多くの会社員は302万程度しかもらえていない。じゃあ、平均年収の400万は何で実現されているかといえば……、年収1000万を超える少数の（おそらく重役たちの）高額所得者に由来しているわけだ。データを見るときには、平均値だけでなく、データの分布、中央値、最頻値などもよく見たほうが良いという好例であろう。