Edited at
RDay 6

統計量の分布をぼーっと眺める 〜中心極限定理観察〜


この記事は何?

 「R Advent Calendar 2018」の第6日目です!

 標本平均、標本中央値、標本範囲中央のヒストグラムを描いて、数学的な話抜きに「眺めて」検討する記事です。

 前半は、Rでヒストグラムを描くにあたってのtipsというか、作業記録的な話になります。

 後半では、母集団(乱数)を取っ替え引っ替えしながら、統計量の分布がどうなるか見ていきます。中心極限定理の可視化、また母集団によっては成り立たないことの可視化でもあります。

 記事のきっかけになったツイートがこちら。


Rでヒストグラム

 ヒストグラムはhist関数で描きます。

 オプションを指定しなかった場合、データの全体が図に入るように描かれます。スタージェスの公式により階級数を決めているそう。

 hist関数のオプションは下のページが参考になります。

 階級の幅を手動調整するには、breaksオプションに境界値を指定してあげます。

# 0から1まで0.1刻み

hist(x, breaks = seq(0, 1, 0.1))


ヒストグラムの一部を描く

 ところで、breaksに指定する境界値はデータの全体をカバーできていないとエラーが出ます。

x <- rnorm(10000)

hist(x, breaks = seq(-4, 4, 0.2))
# hist.default(x, breaks = seq(-4, 4, 0.2)) でエラー:
# 'x' の一部分が数えられていません。多分 'breaks' が 'x' の範囲全体をカバーしていません

 実データを扱うときはいつでもデータの全体を見るべきなのでこの仕様は理解できます。しかし、乱数実験などで分布の中心部や注目部分のヒストグラムを「切り抜いて」描きたいときもあります。

 xlimオプションがそれっぽいと思いきや、

hist(x, breaks = seq(-4, 4, 0.2), xlim = c(-4, 4))

# hist.default(x, breaks = seq(-4, 4, 0.2), xlim = c(-4, 4)) でエラー:
# 'x' の一部分が数えられていません。多分 'breaks' が 'x' の範囲全体をカバーしていません

と相変わらずのエラー……。

 マニュアルによれば「Note that xlim is not used to define the histogram (breaks), but only for plotting (when plot = TRUE).」だそうで、これは表示の制御用で求めているものとは違います。

 望みどおりの結果を得るには、データを加工して範囲外の値を落としてからhist関数に渡すといいようです。

# データを切り抜き

x <- x[-4 < x & x < 4]
hist(x, breaks = seq(-4, 4, 0.2))

 準備完了。


乱数実験

 乱数を標本サイズの個数分生成し、標本平均などを計算、これを繰り返して統計量の分布を得ます。

 標本サイズは、n=5、n=20、n=45とします。中心極限定理が成り立つ状況であれば、nが4倍で分散が4分の1(標準偏差が2分の1)、nが9倍で分散が9分の1(標準偏差が3分の1)になるはず。

 繰り返しは各5万回ずつ行います。手元のMacで数秒で(ストレスなく!)完了するという基準で決めました。

 標本平均(mean)以外の統計量は、標本中央値(median)、標本範囲中央(midrange)です。ちなみに定義は以下のとおり。


  • 標本中央値(nが奇数):小さいほうから(n+1)/2番目の値

  • 標本中央値(nが偶数):小さいほうからn/2番目と(n+1)/2番目の値の中点

  • 標本範囲中央:最小値と最大値の中点

 Rは多くの乱数があらかじめ用意されていて、この手の実験をするのに便利ですね。


分布エントリー

 それでは、母集団(乱数)の分布のエントリーです。パンパカパーン!


  1. 正規分布

     みんな大好き標準正規分布$\text{N}(0,1)$。「標本平均が中心位置を最ももっともらしく推定する」と言われます。


  2. 一様分布

     平均0、分散1に標準化した$\text{U}(-\sqrt{3},\sqrt{3})$。下図で端っこのところだけ柱が低くなっているのは、$\pm\sqrt{3}$が階級の内部のため。


  3. 指数分布

     $\text{Exponential}(1)$。非対称分布なので、平均≠中央値。


  4. ラプラス分布

     平均0、分散1に標準化した$\text{Laplace}(0,1/\sqrt{2})$。正規分布と対照的に、標本中央値が平均の最尤推定量です。


dist1.png

 ここからは中心極限定理の前提を満たしません。

 分散が無限大に発散していて、高頻度で外れ値(outlier)を生み出す曲者たちです。俗にファットテールと呼ばれる分布の中でも極端なものとなります。


  1. コーシー分布

     $\text{Cauchy}(0,1)$。平均も分散もない分布として有名です。(ちなみに自由度1のt分布なのです。)


  2. 自由度2のt分布

     これを母集団とするモデルを考えることは普通しないと思いますが、平均あり・分散なしの分布の代表として取り上げます。


  3. α=1のパレート分布

     $\text{Pareto}(1,1)$。非対称で平均・分散のない分布の代表として。


  4. α=2のパレート分布

     $\text{Pareto}(1,2)$。非対称で平均あり・分散なしの分布の代表として。


dist2.png


結果 〈正規分布〉

左から標本平均、標本中央値、標本範囲中央

dist3.png

 標本平均、標本中央値ともにnが大きくなると分布が母平均の近くに集中するようになりました。標本範囲中央も若干分布が縮小しているように見えます。分布の広がりの順は「標本平均<標本中央値<標本範囲中央」です。

 正規母集団では、標本平均も正規分布に従うので中心極限定理の理想形といえます。

 標本範囲中央は無限区間の分布でnが大きくなると悪化すると予想していましたが、十分外れ値が出にくいということか集中が見られました。(予想外の結果は楽しい!)


結果 〈一様分布〉

左から標本平均、標本中央値、標本範囲中央

dist4.png

 標本平均、標本中央値、標本範囲中央ともに分布が母平均の近くに集中するようになりました。分布の広がりの順は「標本範囲中央<標本平均<標本中央値」です。

 標本中央値の広がりが大きめです。真ん中付近の値が特段出やすくはないことが響いているのかも?

 標本範囲中央は分布が急速に縮小して圧勝ともいえる結果でした。


結果 〈指数分布〉

左から標本平均、標本中央値、標本範囲中央

dist5.png

 標本平均、標本中央値、標本範囲中央どれも非対称の分布になりました。母集団の非対称性が反映されています。

 標本平均と標本中央値では分布の縮小が見られますが、集中先が違います。どうも標本中央値は母中央値の$\ln 2 = 0.693\ldots$の近くに集中しているみたいです。標本範囲中央が徐々に右にずれていくのは上に非有界の分布なので当然でしょうか。

 標本平均は見ての通りガンマ分布で、中心極限定理により正規分布に漸近します。


結果 〈ラプラス分布〉

左から標本平均、標本中央値、標本範囲中央

dist6.png

 標本平均、標本中央値ともに分布が母平均の近くに集中するようになりました。標本範囲中央はnが増えてもあまり変化がないように見えます。

 分布の広がりは最尤推定量の標本中央値が一番小さくなりました。分散が有限の状況でも、標本中央値のほうが標本平均より有効になるケースもあるんですね……。


結果 〈コーシー分布〉

左から標本平均、標本中央値、標本範囲中央

dist7.png

 標本平均はnが大きくなっても変化が見られず、標本中央値は母中央値の近くに集中するようになりました。標本範囲中央はnが大きくなると悪化(分布が拡大)します。

 コーシー分布は和の再生性があり、標本平均もコーシー分布に従うのですが、正規分布の場合と違って分布が縮小しないので標本平均は役に立ちません。コーシー分布の中心位置の推定には標本中央値がよさそうです。


結果 〈自由度2のt分布〉

左から標本平均、標本中央値、標本範囲中央

dist8.png

 標本平均、標本中央値ともに分布が母平均の近くに集中するようになりました。分布の広がりは標本中央値のほうが小さめです。標本範囲中央はコーシー分布のときと同様、nが大きくなると悪化します。

 分散が発散していて中心極限定理の成り立たない状況ですが、標本平均は見た目一致性がありそうな感じです。


結果 〈α=1のパレート分布〉

左から標本平均、標本中央値、標本範囲中央

dist9.png

 標本平均と標本範囲中央は右方向にシフトするような挙動でした。さらにnを増やすと図の範囲からほとんど出て行きそうです。標本中央値は分布が母中央値の2の近くに集中するようになりました。


結果 〈α=2のパレート分布〉

左から標本平均、標本中央値、標本範囲中央

dist10.png

 標本平均は母平均の2の近くに、標本中央値は母中央値の$\sqrt{2} = 1.414\ldots$の近くに集中するようになりました。標本範囲中央はα=1のときと同様に右方向にシフトする感じです。

 中心極限定理の成り立たない状況ですが、標本平均は一致性がありそうに見えます。ただし、分布の広がりは標本中央値より大きめで、縮小も鈍くなっています。


コード

 正規分布のもの。

par(mfrow = c(3, 3))

rep <- 50000
means <- numeric(rep)
medians <- numeric(rep)
midranges <- numeric(rep)

for (i in 1:rep) {
# 正規分布
r <- rnorm(5)
means[i] <- mean(r)
medians[i] <- median(r)
midranges[i] <- (min(r) + max(r)) / 2
}
means <- means[-2 < means & means < 2]
medians <- medians[-2 < medians & medians < 2]
midranges <- midranges[-2 < midranges & midranges < 2]
hist(means, breaks = seq(-2, 2, 0.1), main = "Sample means (n = 5)")
hist(medians, breaks = seq(-2, 2, 0.1), main = "Sample medians (n = 5)")
hist(midranges, breaks = seq(-2, 2, 0.1), main = "Sample midranges (n = 5)")

for (i in 1:rep) {
# 正規分布
r <- rnorm(20)
means[i] <- mean(r)
medians[i] <- median(r)
midranges[i] <- (min(r) + max(r)) / 2
}
means <- means[-2 < means & means < 2]
medians <- medians[-2 < medians & medians < 2]
midranges <- midranges[-2 < midranges & midranges < 2]
hist(means, breaks = seq(-2, 2, 0.1), main = "Sample means (n = 20)")
hist(medians, breaks = seq(-2, 2, 0.1), main = "Sample medians (n = 20)")
hist(midranges, breaks = seq(-2, 2, 0.1), main = "Sample midranges (n = 20)")

for (i in 1:rep) {
# 正規分布
r <- rnorm(45)
means[i] <- mean(r)
medians[i] <- median(r)
midranges[i] <- (min(r) + max(r)) / 2
}
means <- means[-2 < means & means < 2]
medians <- medians[-2 < medians & medians < 2]
midranges <- midranges[-2 < midranges & midranges < 2]
hist(means, breaks = seq(-2, 2, 0.1), main = "Sample means (n = 45)")
hist(medians, breaks = seq(-2, 2, 0.1), main = "Sample medians (n = 45)")
hist(midranges, breaks = seq(-2, 2, 0.1), main = "Sample midranges (n = 45)")

 残りの分布は乱数生成部分のみ。

    # 一様分布

r <- runif(5, -sqrt(3), sqrt(3))
# 指数分布
r <- rexp(5)
# ラプラス分布
r <- (rexp(5) - rexp(5)) / sqrt(2)
# コーシー分布
r <- rcauchy(5)
# 自由度2のt分布
r <- rt(5, df = 2)
# α=1のパレート分布
r <- runif(5)^(-1)
# α=2のパレート分布
r <- runif(5)^(-(1 / 2))


まとめ


  • 標本平均は、平均と分散が存在する(中心極限定理が成り立つ)状況なら、母集団の中心位置の指標として安定した性能を持つ。今回の記事で見た単峰性分布の場合、小さめの標本(nが2桁以下)でも標本平均が正規分布に近づいていくことを確認できた。しかし、混合正規分布のような多峰性分布ではもっと大きな標本が必要かもしれない。

  • 分散が発散している分布で中心極限定理は成り立たない。コーシー分布やα=1のパレート分布のような平均のない分布では、標本平均が母集団の位置を知るのに役に立たないことが分かった。平均あり・分散なしの自由度2のt分布やα=2のパレート分布では、nを大きくすると母平均の近くに集中するのを観察できた。一般的には分からないが、標本平均が無意味とはいえなさそう。

  • 標本中央値にも中心極限定理に類似する法則(母平均ではなく母中央値まわりへの集中)が成り立っている。さらに、中心極限定理の成り立たない状況でもその分布が縮小していくのを観察できた。

 「標本中央値は頑健(robust)だ」とよく言われていますが、今回の実験でもそれが際立つ結果となりました。素性不明の母集団が相手でも意味のある位置の指標となりそうなのは心強いです。


おわりに

 当日に最初書きかけの状態で記事をアップしてしまい、ご迷惑おかけしました。

 もう少し計画的に準備しないといけませんね……。

 もし余裕があれば派生記事を別のアドカレ向けに書こうと思います。(今年、数学カレンダーがないんですよね。Adventarで適当そうなのを探してみよう。。)

 では、引き続き「R Advent Calendar 2018」をお楽しみください!