5 분 소요

Mar 4 (Tue)

교수님이 정말 젊으시다.. 30대 초반의 통계13 출신 교수님! 거의 조교님들이랑 비슷한 연배이실 것 같은데, 벌써 교단에 선지 3년은 넘으신 듯 하다. 참 세상에 천재는 많은가 보다.

이 수업은 R을 통해 여러 수치 모델을 만들고 모의실험(시뮬레이션)을 돌리는 과목이기에, 기본 이상의 R 실력은 알아서 갖춰와야한다고. 파이썬 실력이 있다면 대체는 가능하다시지만, 최대한 R을 열심히 연습해보자. R, Rstudio, Rmarkdown을 활용하여 수업이 진행되고, 과제는 RMD와 워드 파일을 모두 제출해야하는 듯 하다.

칼절평이라고 하시니 놓치는 것 없이 잘 따라가보자. 지난 학기에 비해서 시험의 보너스 점수는 +20에서 +10으로 줄었고, take home 시험도 36-48시간으로 축소하시겠다고. 이래저래 저번 학기보다 학점 따기는 어려워진 것 같기는 하다.

Mar 6 (Thu)

시뮬레이션(모의실험)이 중요한 이유는 무엇일까? 사람 손으로는 전부를 계산하기가 불가능한 문제라든지, 이론이 실제로 현실에서 증명이 가능한지를 보이는데 유용하기 때문이다. 다음 두 개의 예제를 보자. (앞으로 총 4개의 예제를 볼 것이다)

Simulation for Expectation

다음은 수리통계학회(IMS) 회보에 실린 문제이다.

정상적인 주사위를 여러번 굴릴 때, n번 굴렸을 때 나온 눈의 총 합을 Sn이라 하자. Sn이 처음으로 소수(prime number)가 될 때까지 평균적으로 얼마나 주사위를 굴려야하는가? (DasGupta, 2017)

상당히 짧고 단순해보이는 문제이다. 정상적인 주사위라면 항등분포를 따르는 독립시행으로 던져질 것이다. 예를 들어 처음 던지자마자 소수가 될 확률은 1/2이다. 소수 2, 3, 5가 나오면 바로 끝나기 때문이다. 그렇다면 나머지 1, 4, 6이 나오면 최소 한 번은 더 던져야된다. 이런 식의 시나리오를 계산하면 2회차에서 끝날 확률은 2/9(=8/36), 3회차에서 끝날 확률은 5/54(=20/216)이 된다.

DasGupta1

따라서 정답은

\[E(\tau)=1\times\frac{1}{2}+2\times\frac{2}{9}+3\times\frac{5}{54}+\cdots.\]

그런데 이 과정은 무한하다. 계속해서 뻗어나가는 트리 서치(tree search)의 상황을 어떻게 손계산으로 할 수 있을까? 불가능하다. 그러나! 이 강의를 배우는 이유는 시뮬레이션(모의실험)을 통해 이를 해결할 수 있다.

Prime = function(){
  Ind=0; tau=0; S=0;
  while (Ind<1){
    tau=tau+1; r=sample(6,1);
    S=S+r; Ind=is_prime(S)}
  return(tau)}

사용자 지정 함수 Prime은 indicator Ind가 1이 되면(눈의 합이 소수가 되면) 즉시 while문을 빠져나온다. tau는 시행 횟수, r은 매번 랜덤하게 추출한 눈을 가르킨다. 누적합 S에 r을 더한 후 S가 is_prime()에서 소수로 판명되는 순간 지시자가 1이 되며 반복문이 종료된다. (이제 이걸 혼자서 레퍼런스 없이 작성할 수 있어야한다!)

library(primes)
set.seed(3); nsim = 1e4; taus = replicate(nsim,Prime())
### nsim을 점점 늘리다보면 1/2, 2/9, 5/54 및 2.43에 근사해짐
round(table(taus)/nsim,3)
E.taus = cumsum(taus)/(1:nsim)
plot(100:nsim,E.taus[100:nsim],type="l",xlab="Number of experiments",
     ylab="Expectation of tosses",ylim=c(2,2.5),main="")
abline(h=2.43, col="red")

이제 is_prime()을 가져오기 위해 소수 라이브러리를 가져오자. 시뮬레이션은 10^4=1만 번 진행하겠다. 여기서 replicate() 함수는 1만 번 Prime() 함수를 돌린 뒤 그 결과를 배열로써 벡터 taus에 저장하는 역할을 한다.

### replicate 함수 사용법
### n은 시행 횟수 expr는 반복할 함수, simplify=T이면 배열로, F이면 리스트로 저장
replicate(n, expr, simplify = T)

이제 taus에는 만 개의 ‘주사위를 굴린 횟수’ 원소가 들어간 배열이 저장되어있다.

taus
#[1]  1  1  4  3  1  1  4  2  4  3  3  1  1  2  4  1  3  1  1  1  2  2  1  1
#[25]  2  1  1  3  1 (이하 생략)

이 만 개의 데이터를 보기 좋게 ‘주사위를 굴린 횟수’ 별로 상대도수표를 만들어 정리해보자. taus를 table 함수에 통과시켜 분할표의 형태로 정리한 후, nsim(=10000) 횟수로 나누어 상대도수로 처리한다. 이후 소숫점 세 번째 자리에서 반올림하자.

taus
#    1     2     3     4     5     6     7     8     9    10    11    12    13 
#0.502 0.221 0.095 0.062 0.039 0.021 0.016 0.011 0.008 0.007 0.005 0.004 0.003 
#   14    15    16    17    18    19    20    21    22    25    27    29 
#0.002 0.001 0.001 0.001 0.001 0.001 0.000 0.000 0.000 0.000 0.000 0.000

깔끔한 상대도수표가 나왔다. 이번엔 ‘평균 굴린 횟수’를 쭉 줄세워보자. cumsum() 함수를 사용해 매 회차에서 지금까지 기록된 ‘주사위를 굴린 횟수’(=taus)의 평균을 구한 후, E.taus에 저장하자. 1:nsim을 통해 매 횟수별로 분모를 1씩 키워가며 평균을 구할 수 있다.

### cumsum 함수 용례
### 인자로 배열을 넣어주면 각 회차별 누적합을 벡터로 표출한다
cumsum(1:10)
#  [1]  1  3  6 10 15 21 28 36 45 55

이제 100회부터 1만 회까지 E.taus(=주사위를 굴린 평균 횟수)를 쭉 표시해보면 점차 2.43에 수렴해감을 확인할 수 있다. 수식으로 표현하자면

\[\lim_{n \to \infty} \frac{1}{n}\sum_{i=1}^{n}X_i=E[X]\]

DasGupta10k

비교를 위해 nsim=1e5로 바꾸어 10만 회로 실험 횟수를 늘려보자. 여전히 2.43에 수렴함을 볼 수 있다.

DasGupta100k

10만 회를 반복했을 때의 상대분포표를 보면 종전에 우리가 구한대로 1/2, 2/9, 5/54에 더욱 근접해감을 확인할 수 있다. 1만 회때의 상대분포표와 비교해보자. 더 정확해졌다.

taus
#    1     2     3     4     5     6     7     8     9    10    11    12    13 
#0.500 0.222 0.093 0.062 0.038 0.022 0.015 0.012 0.009 0.007 0.005 0.004 0.003 
#   14    15    16    17    18    19    20    21    22    23    24    25    26 
#0.002 0.001 0.001 0.001 0.001 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 
#(이하 생략)

Simulatio for Variance

이번엔 수식으로 구하기 힘든 통계량을 모의실험을 통해 유추해볼 것이다. 다음과 같은 예를 살펴보자.

iid인 포아송분포에서 표본을 추출하여 모수 theta를 추정하고 싶다. 모수 추정은 표본 평균 $$\hat\theta _1=\bar Y$$ 또는 표본 분산 $$\hat\theta _2=S^2$$을 통해 가능하다. 둘 중 무엇이 더 좋은 추정량인가?

여기서 ‘더 좋다’의 개념은 ‘무엇이 더 모수쪽에 집중되어 분포하는가’를 의미한다. 익히 알려져있듯, 포아송분포에서 표본평균과 표본분산은 모두 불편추정량(unbiased estimator)이다.

\[Y\sim Poi(\theta)\rightarrow E(Y)=Var(Y)=\theta\]

즉 둘 다 사용 가능하다. 다만 둘 중 무엇의 성능이 더 좋은지 답하기 위해서는 약간의 수식이 필요하다.

우리는 더 작은 MSE를 갖는 쪽이 더 좋은 추정량이라고 판단할 것이다. MSE는 ‘편향 제곱 더하기 분산’으로 변형이 가능하다. 이를 통해 모수에 관한 관계식이 도출된다.

\[MSE(\hat\theta)=E[(\hat\theta-\theta)]^2==[E(\bar\theta)-\theta]^2+Var(\hat\theta)=Bias(\hat\theta)^2+Var(\hat\theta)\]

그런데 둘 다 불편추정량이라 하였으므로 Bias=0이 된다. 따라서 분산만을 비교하면 된다. 문제는 표본평균의 분산 \(Var(\hat\theta _1)\)은 쉽게 구해지지만, ‘표본분산의 분산’ \(Var(\hat\theta _2)\)은 구하기가 어렵다. 여기서 시뮬레이션이 나온다. 모수 theta를 3이라고 치자.

set.seed(10); theta <- 3 #true value of parameter
poisson.varest <- function(n, theta, iter) {
  theta.1 <- theta.2 <- numeric(iter)
  for(i in 1:iter) {
    Y <- rpois(n, theta); theta.1[i] <- mean(Y); theta.2[i] <- var(Y)}
  print(cbind(var.1=var(theta.1), var.2=var(theta.2)))
  return(cbind(theta1=theta.1, theta2=theta.2)) }
est <- poisson.varest(n=10, theta=theta, iter=50000)

먼저 사용자정의함수 poisson.varest(포아송 분산 추정량)를 정의하자. n은 표본(샘플)의 크기(개수), theta는 모수(참값 3), iter는 반복횟수이다. theta.1과 theta.2에는 길이가 iter인 숫자형 벡터가 저장되고, 0으로 초기화된다. rpois() 함수를 통해 총 5만 번 Poi(3)에서 10개의 샘플을 추출하여 Y에 저장한 후, 각각의 평균과 분산을 벡터 theta.1과 theta.2에 저장한다.

이제 반복문을 나와 theta.1(표본평균)의 분산을 var.1 벡터에, theta.2(표본분산)의 분산을 var.2 벡터에 저장한 후 cbind() 함수를 통해 두 열벡터를 가로 방향으로 묶어서 출력해준다. 또한 이 값들을 묶어 return한 값을 est에 저장한다.

est <- poisson.varest(n=10, theta=theta, iter=50000)
#         var.1    var.2
#[1,] 0.2992606 2.306408

확연히 표본평균의 분산 var.1(=0.299)가 표본분산의 분산 var.2(=2.306)보다 작음을 확인할 수 있다. 더 오밀조밀하게 모인 ‘표본평균’이 그렇지 못한 ‘표본분산’보다 더 좋은 불편추정량임을 알아냈다! 이를 시각화해보면

par(mfrow=c(1,2))
hist(est[,1],freq=FALSE,xlab=expression(hat(theta)[1]),col="gray",main=" ")
hist(est[,2],freq=FALSE,xlab=expression(hat(theta)[2]),col="gray",main=" ")

better_poission

한눈에 보이는 장면에 현혹되지 말자. x축을 보면 표본평균의 분산이 더 좁음을 알 수 있다!

출처: 통계계산방법(STAT323) ㅎㄷㄱ 교수님

참고