[손으로 푸는 통계 ver1.0] 88. 표본분산의 분포 시뮬레이션 (1) 확률밀도함수 비교

 

 

우리가 표본분산의 분포를 유도할 때 설정했던 두가지 조건입니다. 

 

1) 표본평균의 분포가 정규분포를 따를 만큼 표본의 크기 n이 크다. 
2) 모집단의 분포는 정규분포를 따른다. 

 

1번은 표본의 크기를 충분히 크게 하면 되는거구요. 두번째 조건도 표본의 크기가 충분히 크면 무시할 수 있다는 것을 지난시간에 다뤘습니다. 증명하진 않고 증명이 되어 있는 논문만 보여드렸습니다. 

 

오늘은 통계 프로그램인 R을 이용해서 정말 표본의 크기가 충분히 크면 모집단이 정규분포를 따르지 않아도 표본분산이 카이제곱분포를 따르는지 확인해보려고 합니다. 

 

모집단은 1부터 10까지의 자연수로 설정했습니다. 전혀 정규분포가 아닙니다. 

 

모집단 = {1,2,3,4,5,6,7,8,9,10}

 

크기가 2인 표본

크기가 2인 표본을 10000개 뽑아서 각 표본마다 아래 값을 구했습니다. 복원추출을 사용했습니다.

 

$\frac{n-1}{\sigma^{2}}s^{2} $

 

표본의 개수를 10000개 뽑았다는 것은 무한히 많이 뽑는 상황을 가정한 것입니다. 표본의 크기가 아니라 개수입니다. 개수에요. 개수. 개수개수개수개수. 한 표본의 '크기'는 10입니다. 표본 하나를 뽑았을 때 원소가 10개라는 겁니다. 그런 표본을 10000개 뽑은거죠. 

 

히스토그램을 그렸습니다. 전체 넓이가 1이 되도록 히스토그램을 그리는 기능을 사용했습니다. 전체넓이가 1이 되므로 히스토그램은 '확률 밀도 함수' 가 됩니다. 

 

크기가 2인 표본분산을 이용하여 계산된 $\frac{n-1}{\sigma^{2}}s^{2} $ 의 분포는 자유도 1인 카이제곱분포입니다. 좌측에는 표본을 직접 추출해서 만든 히스토그램을 그리고 우측에는 자유도 1인 카이제곱분포 함수를 그렸습니다. 

 

 

표본의 크기가 겨우 2개인데요. 모양이 비슷합니다. 

 

표본의 크기를 10으로 늘려봅시다. 

 

크기가 10인 표본

아래는 크기가 10인 표본의 히스토그램과 9자유도 카이제곱분포 그래프입니다. 빨간 선은 최댓값의 위치입니다. 모양은 비슷한데 최댓값 위치가 다릅니다. 

 

 

랜덤 추출이라서 여러번 더 추출해봤지만 최댓값 위치는 같아지지 않았습니다. 

 

크기를 30으로 늘려봅시다. 

 

크기가 30인 표본

아래 왼쪽 그림은 크기가 30인 표본의 히스토그램입니다. 오른쪽 그림은 29자유도 카이제곱분포입니다. 

 

 

모양은 비슷하고 최댓값 위치도 상당히 비슷해졌습니다. 30이라는 숫자가 참 놀라운 것 같아요. 중심극한정리에서도 표본크기 30이면 정규분포와 얼추 비슷해졌었는데요. 표본분산의 분포에서도 표본크기 30에서 분포가 상당히 비슷합니다. 일반적으로 t검정을 할수 있느냐 없느냐의 기준은 표본 크기 30입니다. 물론 29,28,27 다 비슷하겠지만 30으로 딱 떨어지니 30이라는 숫자를 사용한 것 같아요. 

 

크기를 50으로 키워봅시다. 

 

크기가 50인 표본

아래 왼쪽 그림은 크기가 50인 표본의 히스토그램입니다. 오른쪽 그림은 49자유도 카이제곱분포입니다. 

 

 

최댓값 위치도 일치합니다. 

 

한가지 문제가 모입니다. 전반적인 모양은 일치하지만 함수 값이 다릅니다. 위 그림에서 히스토그램의 최댓값은 0.06 근처인데, 카이제곱분포는 0.04 근처입니다. 최댓값은 히스토그램 구간 간격에 따라서도 변하기 때문에 일단 넘어가겠습니다. 다음시간에는 정량적인 평가를 해보려고 합니다. 위 두 그래프에서 구한 p값을 비교해보겠습니다. 

 

 

<R 코드>

#설정값
ppltn=c(1,2,3,4,5,6,7,8,9,10) 
size=30


#분포함수 그리기
var_p=var(ppltn)*(length(ppltn)-1)/length(ppltn)

dt=list() 

for (i in 1:10000){
dt[[i]]=var(sample(ppltn,size,replace=TRUE))*(size-1)/var_p

}


par(mfrow=c(1,2))

h1=hist(unlist(dt),plot=FALSE)


#x=seq(0,h1$breaks[length(h1$breaks)],by=(h1$breaks[2]-h1$breaks[1]))
#크기 30이상은 아래 x로 설정
x=seq(0,h1$breaks[length(h1$breaks)],by=1)


hist(unlist(dt) ,breaks=x  ,freq=FALSE,xlab="x^2",
     main=paste0( "히스토그램 (표본크기 " , size ,")" ) )
points(c(size-3,size-3),c(0,1),type='l',col='red')

x2=seq(0,h1$breaks[length(h1$breaks)],by=0.01)

plot( x2, dchisq(x2,df=size-1),type='l',ylab='Density',xlab="x^2",
      main=paste0('카이제곱분포 (',size-1,"자유도)" ) )
points(c(size-3,size-3),c(0,1),type='l',col='red')


(max(h1$density)-dchisq(size-3,size-1))/max(h1$density)*100