Order statistics simulation

Summary

본 내용은 STAT-110 강의 25번의 일부 내용을 R 코드로 simulation 합니다. 자세한 내용은 해당 강의를 참조하세요.

\(U1, U2, ,,, Un \sim Unif(0, 1)\) 을 따른다고 하자. 그리고 이 확률변수값을 가장 작은 것에서 가장 큰 것 순서로 정렬을 한다. 그럼, 이중 5번째 위치한 확률변수값의 분포는 어떻게 될까? 즉, 순서 통계량 Order statistics 에 대한 분포를 구한다.

강의 25번에서 이론적인 분포를 아래와 같이 구했다.

\(f_{X(j)} (x) = n {n-1 \choose j-1} F(x)^{j-1} (1 - F(x))^{n-j} f(x)\)

\(Unif(0, 1)\) 의 확률밀도함수 f(x) = 1 이고, CDF 는 선형적으로 증가하는 함수 x 이지.

따라서 위의 공식을 Unif 의 PDF, CDF 로 대체하면 아래와 같다.

\(f_{X(j)} (x) = n {n-1 \choose j-1} x^{j-1} (1 - x)^{n-j}\)

위의 공식은 Beta 분포의 PDF 형태와 같지? \(U(j) \sim Beta(j, n-j+1)\) 이 되는거야.

이 내용을 R 코드로 simulation 하고, Beta 분포를 가시화 해서 확인해보자.

set.seed(1234)

n <- 20 # 20개의 확률변수
j <- 5 # 5번째 순서 통계량 위치

# order(u): 정열된 값에 대한 인덱스를 반환
# order(u)[5]: 정렬된 인덱스의 5번째
# u[order(u)[5]]: 5번째 인덱스의 값 조회
num_trials <- 10000 # simulation 시도 횟수 
trial_result <- sapply(1:num_trials, function(i) {
  u <- runif(n); u[order(u)[j]]
})

# simulation 평균값
mean(trial_result)

## [1] 0.2381767

# simulation 값의 분포를 확인한다.
hist(trial_result, freq = FALSE, breaks = 100, 
     main = paste(j,'th Order statistics simulation', sep=''))


# 이론적인 Beta 분포를 그려서 위의 히스토그램하고 비교한다.
x <- seq(0, 1, by=0.01)

a <- j
b <- n-j+1

# 아래 2 문장은 똑같은 내용임
#y <- (1/beta(a, b))*x^(a-1)*(1-x)^(b-1)
y <- dbeta(x, shape1 = a, shape2 = b) # Beta 분포의 확률밀도함수 사용

# 확률밀도함수 시각화
lines(x, y, type='l', col='red')

# beta 분포 기대값
a/(a+b)

## [1] 0.2380952