順序統計量とは

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

はじめにサイコロの出る目についての問題を考え、そのあとで目標である一様分布の場合を考えます。サイコロの問題は一様分布の離散版と考えることができます。

サイコロの場合

4個のサイコロを振って小さい順から2番目の目が2以下である確率を計算してみましょう。これは、4個全てが2以下、4個中3個だけが2以下、4個中2個だけが2以下という3つの排反なパターンに分けることができます。それぞれの確率は

(2/6)⁴
4 * (2/6)³ * (4/6)
4C2 * (2/6)² * (4/6)²

となり、合計すると11/27となります。

同様に4個中2番目がx以下である確率は

(x/6)⁴
4 * (x/6)³ * (1-x/6)
4C2 * (x/6)² * (1-x/6)²

となり、合計するとx²(3x²-8x+6)/6⁴となります。

さらに一般化して、4個中p番目がx以下である確率は、4個中4個がx以下から4個中p個だけがx以下までの4-p+1通りの排反なパターンに分けることができます。4個中k個だけがx以下である確率は

4Ck * (x/6)^k * (1-x/6)^4-k

となり、合計すると $\displaystyle{\sum_{k=p}^4 {}_4C_k \frac{x}{6}^k (1-\frac{x}{6})^{4-k}}$ となります。

最後にn個中p番目がx以下の場合もやっておきましょう。n個中p番目がx以下である確率は、n個中n個がx以下からn個中p個だけがx以下までのn-p+1通りの排反なパターンに分けることができます。n個中k個だけがx以下である確率は

nCk * (x/6)^k * (1-x/6)^n-k

となり、合計すると

$$ \sum_{k=p}^n {}_nC_k \frac{x}{6}^k (1-\frac{x}{6})^{n-k} $$

となります。

一様分布の場合の例

0以上1以下の一様分布に従う独立な確率変数X_1, X_2, X_3, X_4の2番目の値が0.3以下である確率は、

4個全て0.3以下であるのが0.3⁴
4個中3個だけが0.3以下であるのが4 * 0.3³ * 0.7
4個中2個だけが0.3以下であるのが4C2 * 0.3² * 0.7²

であり、それの合計です。

同様に4個中2番目がx以下である確率は、

4個全てx以下であるのがx⁴
4個中3個だけがx以下であるのが4x³ (1-x)
4個中2個だけがx以下であるのが4C2 x² (1-x)²

となり、これの和が分布関数です。

確率密度関数はこれのxでの微分です。和をとる前に各項の微分を計算しておくと、それぞれ

4 x³
4 3x² (1-x) + 4 x³ (-1) = 4 3x² (1-x) - 4 x³
4C2 * 2x (1-x)² + 4C2 x² * (-2)(1-x) = 4C2 * 2 x (1-x)² - 4 * 3 x² (1-x)

となります。これを足すと、打ち消しあう項があるので、

$$ 4 x^3 + (4 \times 3x^2 (1-x) - 4 x^3) + ({}_4C_2 2x(1-x)^2 - 4 \times 3 x^2 (1-x)) = {}_4C_2 2x(1-x)^2 $$

となります。

この結果をシミュレーションと比較してみましょう。

library(dplyr)
library(ggplot2)

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 <- 4
p <- 2
fun <- function(x) 12*x*(1-x)**2
compare(N,n,p,fun)

を実行すると、

f:id:unaoya:20190104204457p:plain — n=4, p=2の密度関数との比較

となります。

一般的な状況

最後により一般に、一様分布に独立に従う標本の順序統計量の確率密度関数を計算してみます。

0以上1以下の一様分布に従う独立な確率変数X_1, ..., X_nのp番目の値がx以下である確率は、

n個全てx以下であるのがxⁿ
n個中n-1個だけがx以下であるのがn x^n-1 (1-x)
n個中k個だけがx以下であるのがnCk x^k (1-x)^n-k
n個中p個だけがx以下であるのがnCp x^p (1-x)^n-p

の和が分布関数です。

上と同様に確率密度関数はこれのxでの微分です。和をとる前に各項の微分を計算しておくと、それぞれ

n x^n-1
n(n-1) x^n-2 (1-x) - n x^n-1
nCk k x^k-1 (1-x)^n-k - nCk(n-k) x^k (1-x)^n-k-1
nCp p x^p-1 (1-x)^n-p - nCp(n-p) x^p (1-x)^n-p-1

となります。これを足すと、打ち消しあう項があるので、

$$ n x^{n-1} + (n(n-1)x^{n-2} (1-x) - n x^{n-1}) \\ + ({}_nC_2 (n-2)x^{n-3}(1-x)^2 - n(n-1)x^{n-2} (1-x)) \\ + \cdots \\ + ({}_nC_p px^{p-1}(1-x)^{n-p} - {}_nC_p(n-p)x^p (1-x)^{n-p-1}) \\ = {}_nC_p px^{p-1}(1-x)^{n-p} $$

となります。これがいわゆるベータ分布です。

前回冒頭で紹介したケースでシミュレーションと比較してみましょう。サンプルサイズn=10でp=3とすると、密度関数は上の計算結果から

$$ {}_{10}C_33x^2 (1-x)^7 $$

となります。

library(dplyr)
library(ggplot2)

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 <- 10
p <- 3
fun <- function(x) (10*9*8)/2*x**2*(1-x)**7
compare(N,n,p,fun)