筆者の学生時代¹、数学で統計の分野を触ったときに、このようなことを教わりました。

母集団からいくつかの標本をピックして、そこから母集団の傾向を推測する (標本調査) とき、

確率的に標本の平均は母集団の平均をよく表すが、

標準偏差のような**『ばらつき』を示す尺度**は、標本のそれは母集団のそれよりも確率的に少なくなりがちである。

要するに、

標本の標準偏差は母集団より小さくなるんだよ! ΩΩΩ

ということです。

その後すぐに標本数 $n$ ではなく $n - 1$ で割る**「不偏分散」**の式を教わったりしたものですが、まあ当時は「ほーん」と聞き流していたその情報、確かめたいと思ったことはないでしょうか?

ちゃんとした代数学の証明で「不偏分散の期待値は母集団の分散に等しい」なーんて言われても、確かにそれはそうなのでしょうが、やはり具体的なデータ²を見ないといまいちスッキリしない、なんて経験はないでしょうか?

とはいっても、それを確かめるためのデータなんてありませんし、標本の偏りを考えれば一度や二度計算しただけでは確信など持てるはずもなく……。

ちょっと待ってください。あなたの前にある箱 (あるいは板) はそのためにあるのではないでしょうか? ~~違う気もします。でも少なくとも 5ch で煽り合いをするためにあるものではないはずです。~~ コンピュータがあればコインを一万回投げるのだって朝飯前です。

子供時代にしまい込んだ素朴な疑問、計算機科学で殴りつけてやりましょう。

仮想世界でイキる

今回、シミュレーションの実行には Python + JupyterLab を使用します。が、ぶっちゃけなにを使おうが構わないでしょう。最近私が妙に Python の記事を Qiita に投稿しつづけているから以上の意味はありません。だからプログラミング言語でマウント合戦とか頭の悪いことは今すぐやめろ。自然言語ならなおさらです³。

In[]

import math

import matplotlib.pyplot as plt
import numpy as np

今回母集団は平均 50、標準偏差 10 の正規分布に従う集団を想定します。余談も余談なのですが、この $\mathcal{N}(50, 10^2)$ は偏差値の定義でもあります。一般には「平均が 50 になるように変形した」という部分しか知られていませんが、標準偏差 10 という部分がかなり重要です。これがあるおかげで「その成績がどの程度抜きん出ているのか」が統一的に議論できるわけですから。おまけに 50 なんて意味ありげな数字選んでるせいで下限が 0 で上限が 100 なんてとんでもない勘違いしている人も多い始末⁴。

In[]

mu = 50.
sigma = 10.

抽出する標本は 10 個とします。母集団の規模にもよりますが、これは調査としてはかなり少ないと思います。トリビアの種でも最低 2000 人とか言われていたこのご時おい今回の記事脱線多すぎと違うか? これだけ少ない標本でも期待値で見たときにちゃんと母集団と一致しているのか楽しみです。

In[]

samples_num = 10

平均

まずは「標本調査でも (確率的には) そのまま使っても問題ない」と言われている指標である、平均を見てみましょう。

In[]

np.random.seed(555)

mean_list = []
for _ in range(10000):
    samples = np.random.normal(mu, sigma, samples_num)
    mean_list.append(np.mean(samples))
mean = np.array(mean_list)

plt.hist(mean, bins=30, density=True, histtype='stepfilled')
plt.title(F'Mean: $\\mu = {np.mean(mean):.3f}$, $\\sigma = {np.std(mean):.3f}$')
plt.show()

結果、10000 回の標本調査の結果の平均は 50.022 となりました。母集団の平均が 50 なので、ちゃんと推定できています。標準偏差 3.218 は多いのか少ないのかはわかりませんが、おそらくは標本を増やせばばらつきは小さくなっていくのでしょう。

標準偏差

つづいて本日の主役、標準偏差くんです。

In[]

np.random.seed(913)

std_list = []
for _ in range(10000):
    samples = np.random.normal(mu, sigma, samples_num)
    std_list.append(np.std(samples))
std = np.array(std_list)

plt.hist(std, bins=30, density=True, histtype='stepfilled')
plt.title(F'Standard Deviation: $\\mu = {np.mean(std):.3f}$, $\\sigma = {np.std(std):.3f}$')
plt.show()

結果、10000 回の標本調査の平均は 9.235 となりました。母集団の標準偏差が 10 なので、確かに少なく見積もられています。

標準偏差の不偏推定量⁵

標準偏差とは分散の (非負な) 平方根です。分散とは各データの**「平均との差を二乗した値」の平均**です。
$$
s = \sqrt{s^2} = \sqrt{\frac{1}{n}\sum_{i = 1}^n(x_i - \overline{x})^2}
$$
標本を無作為に選ぶときに、平均からあまりにも外れた大きい値や小さい値が選ばれる可能性は低いですから (一様分布みたいなばらつきもクソもねえ分布はともかく)、どうしたって低く見積もられてしまいます。

なんでそうなるかは各種統計についてもっと真面目に書かれた書物やウェブサイトを参照してもらうとして、$n$ で割って平均を求めるのではなく $n - 1$ で割って得られる値は、その期待値が母集団の分散に一致します。
$$
u^2 = \frac{1}{n - 1}\sum_{i = 1}^{n}(x_i - \overline{x})^2 = \frac{n}{n - 1}s^2
$$

$$
E(u^2) = \sigma^2
$$

じゃあその平方根の期待値は母集団の標準偏差に一致するのかというと微妙に違って、以下の式で求められるそうです。
$$
D = \sqrt{\frac{n - 1}{2}}\frac{\Gamma\left(\frac{n - 1}{2}\right)}{\Gamma\left(\frac{n}{2}\right)}u = \sqrt{\frac{n}{2}}\frac{\Gamma\left(\frac{n - 1}{2}\right)}{\Gamma\left(\frac{n}{2}\right)}s
$$

$$
E(D) = \sigma
$$

なんじゃこの式わ。ちなみに自然数の範囲において $\Gamma(n) = (n - 1)!$ らしいです。 知らねえよ。

ともあれ、この記述を参考に $s$ の係数を求める関数を書いていきます。

In[]

def ustd_coefficient(n):
    try:
        return math.sqrt(n / 2) * math.gamma((n - 1) / 2) / math.gamma(n / 2)
    except OverflowError:
        return math.sqrt(n / (n - 1.5))

In[]

np.random.seed(333)

D_list = []
for _ in range(10000):
    samples = np.random.normal(mu, sigma, samples_num)
    D_list.append(np.std(samples) * ustd_coefficient(len(samples)))
D = np.array(D_list)

plt.hist(D, bins=30, density=True, histtype='stepfilled')
plt.title(F'Unbiased Standard Deviation: $\\mu = {np.mean(D):.3f}$, $\\sigma = {np.std(D):.3f}$')
plt.show()

結果、10000 回の標本調査の平均は 10.024 となりました。母集団の標準偏差が 10 なので、確かに母集団の標準偏差を推定できているようですね。

中央値

ここから先はおまけです。

ちょっとだけ統計をかじった人は**「平均はクソ! 平均依存症から抜け出そう!」みたいな思想に染まるらしいですが^[要出典]、そういう人が決まって持ち出すのが中央値**です。この際なのでこっちもやってみましょう。

正規分布においては平均と中央値は一致するので、母集団の中央値も 50 となります。

私なりにどうなるか予想してみます。まず中央値とは全部の値を 1:1 に分割する値です。で、標本を取り出すときに母集団の中央値を下回る値になる確率と上回る値になる確率は、定義から明らかに等しいとわかります。ですからその標本を 1:1 に分割する値はやはり母集団のそれに収束するでしょう。

では、確かめてみます。

In[]

np.random.seed(753)

median_list = []
for _ in range(10000):
    samples = np.random.normal(mu, sigma, samples_num)
    median_list.append(np.median(samples))
median = np.array(median_list)

plt.hist(median, bins=30, density=True, histtype='stepfilled')
plt.title(F'Median: $\\mu = {np.mean(median):.3f}$, $\\sigma = {np.std(median):.3f}$')
plt.show()

結果、10000 回の標本調査の平均は 50.040 となりました。母集団の中央値が 50 なので、やはり予想どおり、標本から正しく推定できています。

正規四分位範囲

最後に正規四分位範囲をみてみます。前回の記事で確かめたとおり、正規分布においては標準偏差と正規四分位範囲は一致するので、母集団の正規四分位範囲も 10 となります。

あ、NumPy には正規四分位範囲を計算する関数はないので、でっち上げておきます。

In[]

def np_iqr(a, axis=None, out=None, overwrite_input=False, keepdims=False):
    q1 = np.quantile(a, q=0.25, axis=axis, overwrite_input=overwrite_input, keepdims=keepdims)
    q3 = np.quantile(a, q=0.75, axis=axis, overwrite_input=overwrite_input, keepdims=keepdims)
    return np.subtract(q3, q1, out=out)

def np_niqr(a, axis=None, out=None, overwrite_input=False, keepdims=False):
    return np.divide(np.iqr(a, axis=axis, overwrite_input=overwrite_input, keepdims=keepdims), 1.3489795003921634, out=out)

np.iqr = np_iqr
np.niqr = np_niqr

私なりにどうなるか予想してみます。まず、正規四分位範囲は上と下のクォンタイルの差です。クォンタイルの中でも極端なものとして、最小値と最大値を考えてみたとき、標本の最小値は大抵の場合母集団の最小値より大きいでしょうし、標本の最大値は大抵の場合母集団の最大値より小さいでしょう。つまり、標本のクォンタイルは母集団のそれよりも内側に偏りがちになります (中央値はちょうど真ん中のクォンタイルなので影響を受けませんでした)。結果、正規四分位範囲は母集団のそれよりも低く見積もられると予想できます。

では、確かめてみます。

In[]

np.random.seed(315)

niqr_list = []
for _ in range(10000):
    samples = np.random.normal(mu, sigma, samples_num)
    niqr_list.append(np.niqr(samples))
niqr = np.array(niqr_list)

plt.hist(niqr, bins=30, density=True, histtype='stepfilled')
plt.title(F'Normalized Interquartile Range: $\\mu = {np.mean(niqr):.3f}$, $\\sigma = {np.std(niqr):.3f}$')
plt.show()

結果、10000 回の標本調査の平均は 8.640 となりました。母集団の正規四分位範囲が 10 なので、やはり予想どおり少なく見積もられています。

あとがき

正規四分位範囲の不偏推定量ってどうやって求めるんだろう。

専攻は化け学でした。 ↩
数学は具体(実例)より抽象(一般化)に重きをおく学問ですからね、仕方ないね。 ↩
なんで鬼の首とったように英語の不規則性を dis る連中が一定数湧くんでしょうね。気持ち悪いことこの上ないです。 ↩
かくいう私も子供のころは「平均以下を 0 から 50 に変形して、平均以上を 50 から 100 に変形したらすごいいびつな分布になるんじゃ」とかいらぬ心配をしていた一人です。 ↩
このようなやたらと長ったらしい言い回しになっているのはこの分野において用語が混乱しているからです。標本の分散 (標準偏差) を「標本分散 (標本標準偏差)」と呼ぶのか、標本調査に使われるので不偏推定量のほうを「標本分散(標本標準偏差)」と呼ぶのか。不偏標準偏差とは「母集団の標準偏差の不偏推定量」なのか「母集団の分散の不偏推定量の平方根」なのか。あーもうめちゃくちゃだよ。 ↩