標準偏差の不偏推定量はホントに不偏なのか、「コインを一万回投げて」確かめてみた

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

>  母集団からいくつかの標本をピックして、そこから母集団の傾向を推測する **(標本調査)** とき、
>
> * 確率的に**標本の平均は母集団の平均をよく表す**が、
> * 標準偏差のような**『ばらつき』を示す尺度**は、**標本のそれは母集団のそれよりも確率的に少なくなりがち**である。

要するに、

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

ということです。

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

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

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

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

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

[^major]: 専攻は化け学でした。
[^abstraction]: 数学は具体(実例)より抽象(一般化)に重きをおく学問ですからね、仕方ないね。

## 仮想世界でイキる

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

```py:In[]
import math

import matplotlib.pyplot as plt
import numpy as np
```

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

```py:In[]
mu = 50.
sigma = 10.
```

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

```py:In[]
samples_num = 10
```

[^English]: なんで鬼の首とったように英語の不規則性を dis る連中が一定数湧くんでしょうね。 気持ち悪いことこの上ないです。
[^deviation]: かくいう私も子供のころは「平均以下を 0 から 50 に変形して、平均以上を 50 から 100 に変形したらすごいいびつな分布になるんじゃ」とかいらぬ心配をしていた一人です。

### 平均

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

```py: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()
```

![mean.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/241027/0ff32724-a291-ceff-cd67-4b67e1c0f67e.png)

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

### 標準偏差

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

```py: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()
```

![std.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/241027/a07ca591-6f33-1a7e-39cb-ee19205ac215.png)

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

### 標準偏差の不偏推定量[^unbiased]

**標準偏差**とは**分散の** (非負な) **平方根**です。 **分散**とは各データの**「平均との差を二乗した値」の平均**です。
$$
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$ の係数を求める関数を書いていきます。

```py: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))
```

```py: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()
```

![unbiased_std.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/241027/ba72f713-66da-00f1-6ab6-d4395cfcd845.png)

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

[^unbiased]: このようなやたらと長ったらしい言い回しになっているのはこの分野において用語が混乱しているからです。 標本の分散 (標準偏差) を「標本分散 (標本標準偏差)」と呼ぶのか、標本調査に使われるので不偏推定量のほうを「標本分散(標本標準偏差)」と呼ぶのか。 不偏標準偏差とは「母集団の標準偏差の不偏推定量」なのか「母集団の分散の不偏推定量の平方根」なのか。 あーもうめちゃくちゃだよ。

### 中央値

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

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

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

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

では、確かめてみます。

```py: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()
```

![median.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/241027/c063d168-2434-a117-d17f-3ba9357d98f8.png)

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

### 正規四分位範囲

最後に**正規四分位範囲**をみてみます。 [前回の記事](https://qiita.com/BlueRayi/items/ce0f9b3cbb57bfed855f)で確かめたとおり、**正規分布においては標準偏差と正規四分位範囲は一致する**ので、**母集団の正規四分位範囲も 10** となります。

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

```py: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
```

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

では、確かめてみます。

```py: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()
```

![niqr.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/241027/0e99bd37-5c89-2380-0c46-c5bc5e62d885.png)

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

## あとがき

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