ぎるばーとのノート

もっともっと遠くへ行きたい　空が広く見える場所まで

一様分布の順序統計量【ときどき分布・1記事目】

数学・統計学

　ゆるい連載シリーズ「ときどき分布」を始めます。
　シリーズ名のとおり、ときどき確率分布の話をしていこうと思います。

  • 定理の証明は省略
  • 関係式は結果のみ示す
  • 確率密度関数を図示
  • とりあえず乱数実験で確認

といったゆるい感じのスタンスでいきます。

　確率分布の役立つ&役立たない豆知識を紹介していく予定です。
　どうぞお楽しみに！

一様分布の順序統計量

　初回のお題は一様分布の順序統計量です。

一様分布

　一様分布は、分布の台の区間で確率密度関数が一定値である分布です。

X \sim \text{Uniform}(a,\ b)

f(x) = \begin{cases} \frac{1}{b - a} & \text{for $a \leq x \leq b$,}\\ 0 & \text{otherwise.} \end{cases}

　特に、[0, 1]を台とする一様分布を標準一様分布といいます。

X \sim \text{Uniform}(0,\ 1)

f(x) = \begin{cases} 1 & \text{for $0 \leq x \leq 1$,}\\ 0 & \text{otherwise.} \end{cases}

標準一様分布

順序統計量

　標準一様分布からサンプリングした標本（サイズn）について、i番目に小さい観測値をx(i)と名付けます。ここで、x(1)は標本最小値、x(n)は標本最大値です。
　x(1) ... x(n)は、サンプリングを行うたびに異なる値をとります。いったいどんな分布に従うでしょうか？

　観測値x(i)に対応する統計量（順序統計量）をX(i)で表します。X(i)は次の分布に従うことが知られています。

X_{(i)} \sim \text{Beta}(i,\ n + 1 - i)

　たとえば、n = 6であれば、次のようになります。

\begin{cases} X_{(1)} = X_\text{min} &\sim \text{Beta}(1,\ 6) \\ X_{(2)} &\sim \text{Beta}(2,\ 5) \\ X_{(3)} &\sim \text{Beta}(3,\ 4) \\ X_{(4)} &\sim \text{Beta}(4,\ 3) \\ X_{(5)} &\sim \text{Beta}(5,\ 2) \\ X_{(6)} = X_\text{max} &\sim \text{Beta}(6,\ 1) \end{cases}

一様分布の順序統計量（n = 6）

乱数実験

　前節の順序統計量の理論分布は、数学的に導かれたものです。ここでは数学的議論に立ち入ることなく、乱数を使った実験で確認してみます。

  1. 以下を1万回繰り返す（回数に特に根拠はなし）
    1. 標準一様乱数を6個生成
    2. 小さい順に整列
    3. 順序統計量の値をプール
  2. 乱数実験の結果（ヒストグラム）と理論分布を重ねて図示

　乱数実験はRStudioで行いました。Rのコードは記事の末尾（続きを読む以下）にあります。

乱数実験結果

　よく一致する結果が得られました！

付録：Rコード集

乱数実験本体
rep <- 10000
x_1 <- numeric(rep)
x_2 <- numeric(rep)
x_3 <- numeric(rep)
x_4 <- numeric(rep)
x_5 <- numeric(rep)
x_6 <- numeric(rep)

for (i in 1:rep) {
	x <- sort(runif(n = 6))
	x_1[i] <- x[1]
	x_2[i] <- x[2]
	x_3[i] <- x[3]
	x_4[i] <- x[4]
	x_5[i] <- x[5]
	x_6[i] <- x[6]
}
結果図示
library(RColorBrewer)
cols <- brewer.pal(8, "Set2")

par(mfrow=c(2, 3))

hist(x_1, breaks = seq(0, 1, 0.05), probability = T)
curve(dbeta(x, 1, 6), add = T, col = cols[1], lwd = 2)

hist(x_2, breaks = seq(0, 1, 0.05), probability = T)
curve(dbeta(x, 2, 5), add = T, col = cols[2], lwd = 2)

hist(x_3, breaks = seq(0, 1, 0.05), probability = T)
curve(dbeta(x, 3, 4), add = T, col = cols[3], lwd = 2)

hist(x_4, breaks = seq(0, 1, 0.05), probability = T)
curve(dbeta(x, 4, 3), add = T, col = cols[4], lwd = 2)

hist(x_5, breaks = seq(0, 1, 0.05), probability = T)
curve(dbeta(x, 5, 2), add = T, col = cols[5], lwd = 2)

hist(x_6, breaks = seq(0, 1, 0.05), probability = T)
curve(dbeta(x, 6, 1), add = T, col = cols[6], lwd = 2)