の続きで、標本の中央値の分布がどのようになるかを調べていきます。また、より一般に順序統計量というものについて調べます。

順序統計量というのは、標本におけるある順位の値のことです。例えば、最大値や最小値、中央値などがそれに当たります。より一般には、サンプルサイズ10のサンプルのうち3番目のデータがどのような分布になるか？などを考えます。

特に今回と次回は0以上1以下の一様分布から独立なサイズnのサンプルを取った時、そのp番目のデータがどのような分布を持つかを調べていきます。

シミュレーション

まずは実験してみましょう。

例えば、サンプルサイズn=10の標本のうち小さい方からp=3番目の値がどのような分布になるかをシミュレーションしてみます。

set.seed(123)
runif(10)

を実行すると

 [1] 0.2875775 0.7883051 0.4089769 0.8830174 0.9404673
 [6] 0.0455565 0.5281055 0.8924190 0.5514350 0.4566147

と乱数が生成され、これの3番目の値は0.37053701です。

もう一度

runif(10)

を実行すると

 [1] 0.95683335 0.45333416 0.67757064 0.57263340
 [5] 0.10292468 0.89982497 0.24608773 0.04205953
 [9] 0.32792072 0.95450365

と乱数が生成され、これの3番目の値は0.24608773です。

この3番目の値をとる、というのをN=100000回繰り返してヒストグラムを描いてみます。

ord_sampling <- function(N,n,p){
  x <- c()
  for(i in 1:N){
    x[i] <- sort(runif(n))[p]
  }
  data.frame(x=x) %>%
    ggplot(aes(x=x)) +
    geom_histogram(binwidth = 0.05,
                   boundary = 0,
                   fill = 'royal blue') +
    labs(title = paste("n =",as.character(n),", p = ",as.character(p)))
}

runif(10)
N <- 1e+5
n <- 10
p <- 3
ord_sampling(N,n,p)

f:id:unaoya:20181229211425p:plain — 標本3番目の分布

おおよそ0.2から0.25あたりの値が最頻で、直感にはまあ合う結果でしょう。

他にも、同じようにサンプルサイズn=10の標本のうち小さい方からp=1,2,4,5番目の値がどのような分布になるかをシミュレーションしてみます。

f:id:unaoya:20181229211536p:plain — 標本最小値の分布

f:id:unaoya:20181229211539p:plain — 標本2番目の分布

f:id:unaoya:20181229211604p:plain — 標本4番目の分布

f:id:unaoya:20181229211617p:plain — 標本5番目の分布

先に答えを言ってしまうと、実はこれはベータ分布と呼ばれる分布になっています。このことを確率密度関数を計算することで確かめていきましょう。

今回は小手調べでサンプルサイズn=2の時に、最小値の分布がどのような確率密度関数を持つのかを計算してみます。

確率密度関数を計算するには、分布関数を求め、その微分を計算する必要があります。ここでは分布関数は、xに対してサンプルの最小値がx以下である確率を与える関数です。今回は、サンプルX1, X2が一様分布U(0,1)に独立に従う確率変数で考えます。そのうち小さい方（正確には大きくない方）がx以下である確率を計算してみましょう。

サイコロの場合

一様分布の場合を計算する前に、より単純な設定でサイコロの出目についての問題を考えてみます。サイコロ二つを振って小さい方（正確には大きくない方）が2以下である確率を計算してみましょう。

サイコロ二つの名前を1番と2番とし、1番の目をX1とし、2番の目をX2とします。これらは確率変数で、1から6までの値を1/6ずつの確率でとります。 X1が2以下になる確率は2/6で、X2が2以下になる確率も2/6です。

最小値が2以下になる確率は、

X1が2以下でX2が3以上になる確率2/6 * (1-2/6) = 8/36
X1が3以上でX2が2以下になる確率(1-2/6) * 2/6 = 8/36
X1もX2も両方2以下である確率2/6 * 2/6 = 4/36

の合計なので、20/36となります。これは二つの出目の36通りのうち、条件に該当するものがどちらかが1または2である20通りであることからも計算できます。

同様に小さい方がx以下である確率は次のように計算できます。 X1が2以下になる確率はx/6で、X2が2以下になる確率もx/6です。最小値がx以下になる確率は、

X1がx以下でX2がx+1以上になる確率x/6 * (1-x/6) = x/6 - x²/36
X1がx+1以上でX2がx以下になる確率(1-x/6) * x/6 = x/6 - x²/36
X1もX2も両方x以下である確率x/6 * x/6 = x²/36

の合計なので、2x/6 - x²/36となります。この結果はどちらかがx以下であるパターンの数え上げで2xから共通部分のx²を引いたものとしても計算できます。

一様分布の場合

さて、上のサイコロと同様の計算を行います。サンプルとして与えられる二つの値を確率変数X1, X2とします。これらは独立に一様分布U(0,1)に従う確率変数であり、X1がx以下である確率はxであり、X2がx以下である確率も同じくxです。

最小値がx以下になる確率は、

X1がx以下でX2がx以上になる確率x * (1-x) = x - x²
X1がx以上でX2がx以下になる確率(1-x) * x = x - x²
X1もX2も両方x以下である確率x * x = x²

の合計なので、2x - x²となります。これが分布関数です。さらにこれを微分すると2-2xで、これが密度関数です。

最後にシミュレーション結果と、上で求めた密度関数のグラフを重ねて表示してみましょう。

compare <- function(N,n,p,fun){
  x <- c()
  for(i in 1:N){
    x[i] <- sort(runif(n))[p]
  }
  data.frame(x=x) %>%
    ggplot(aes(x=x)) +
    geom_histogram(aes(y=..density..),
                   fill = 'royal blue',
                   binwidth = 0.05,
                   boundary = 0) +
    stat_function(fun=fun) +
    labs(title = paste("n =",as.character(n),", p = ",as.character(p)))
}

N <- 1e+5
n <- 2
p <- 1
fun <- function(x) 2*(1-x)
compare(N,n,p,fun)