[손으로 푸는 통계 ver1.0] 91. 표본분산의 분포 시뮬레이션 (4) p값 비교

 

 

표본분산의 분포와 카이제곱분포를 비교하고 있습니다.

 

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

 

지난시간에는 누적분포함수를 비교했습니다. 모집단이 균등분포를 따르는 경우 표본을 아무리 크게 해도 표본분산의 분포와 카이제곱분포가 일치하지 않았습니다. 

 

이번 글에서는 p값을 비교해봅시다. 정량적인 비교입니다. 

비교 방법을 설명하겠습니다. 카이제곱분포의 좌측꼬리 p값이 0.05가 나오는 확률변수 값을 히스토그램에 적용하여 p값을 구합니다. 0.05와 구해진 p값을 비교하면 됩니다. 우측꼬리에서도 같은 방법으로 p값을 구합니다. 

 

모집단의 종류, 표본의 크기를 바꿔가며 구했습니다. 히스토그램은 표본 10000개를 뽑아서 그렸습니다. 히스토그램을 그리고 해당 히스토그램에서 p값을 구했는데요. 랜덤추출이라서 히스토그램이 조금씩 다르게 그려질 수 있으므로 10회 반복하여 평균과 표준편차를 구했습니다. 표준편차는 0이 나왔습니다. 

 

결과는 아래와 같습니다. (모집단이 1:10 은 1부터 10까지의 자연수를 의미합니다. )

 

 

균등분포를 따르는 모집단의 경우 p값이 카이제곱분포와 많이 다릅니다. 표본의 크기가 커지면 차이가 더 커집니다. 정규분포를 따르는 모집단은 p값이 카이제곱분포와 비슷한 것을 알 수 있습니다. 

 

다음시간에는 지난번에 그냥 넘어갔던 증명을 시도해보겠습니다. 아래 증명입니다. 

 

"모집단의 분포와 상관 없이, n이 무한대로 가면 표본분산의 분포는 카이제곱분포에 가까워져 간다."

 

 

 

<사용 코드>

library(dplyr)

#1.모집단 설정
#ppltn=c(1,2,3,4,5,6,7,8,9,10) 
#ppltn=1:1000
#ppltn=rnorm(10) 
ppltn=rnorm(1000) 

var_p=var(ppltn)*(length(ppltn)-1)/length(ppltn) #모분산

#2. 표본 크기 설정
size=c(30,50,100,500,1000,3000)

#3. 비교 p값 설정(우측꼬리기준으로)
p=0.05


#4. 실험

#set.seed(9999)
df_result=data.frame(n=NA,hist_L_p5=NA,dist_L_p5=NA,diff1=NA,hist_R_p5=NA,dist_R_p5=NA,diff2=NA)

for (k in 1:length(size)){
#for (k in 6){
  
  dt=list() 
  df=data.frame(n=NA,hist_L_p5=NA,dist_L_p5=NA,diff1=NA,hist_R_p5=NA,dist_R_p5=NA,diff2=NA)
  
  
  
  for (j in 1:10){ #30회 반복
    
    for (i in 1:50000){ 
      dt[[i]]=var(sample(ppltn,size[k],replace=TRUE))*(size[k]-1)/var_p 
    }
    
    dt=unlist(dt)
    dt_sort=sort(dt)
    chisq_left=qchisq(p,size[k]-1)
    chisq_right=qchisq(1-p,size[k]-1)
    
    
    hist_L_p5=order(abs(dt_sort-chisq_left))[1]/50000 %>% round(3)
    dist_L_p5=p
    diff1=(hist_L_p5-dist_L_p5)/dist_L_p5*100 
    diff1=diff1 %>% round(2)
    
    hist_R_p5=1-order(abs(dt_sort-chisq_right))[1]/50000 %>% round(3) 
    dist_R_p5=p  %>% round(2)
    diff2=(hist_R_p5-dist_R_p5)/dist_R_p5*100  %>% round(2)
    diff2=diff2 %>% round(2)
    
    
    df[j,]=list(size[k],hist_L_p5,dist_L_p5,diff1,hist_R_p5,dist_R_p5,diff2)
    
  }
  
  
  mean=apply(df,2,mean) %>% round(3)
  sd=apply(df,2,sd) %>% round(2)
  
  
  
  df_result[k,]=list(mean[1],
                     paste0(mean[2]," ± ",sd[2]),paste0(mean[3]," ± ",sd[3]),paste0(round(mean[4],1)," ± ",sd[4]),
                     paste0(mean[5]," ± ",sd[5]),paste0(mean[6]," ± ",sd[6]),paste0(round(mean[7],1)," ± ",sd[7]))
  
}

View(df_result)


write.csv(df_result,file="ex1.csv")

########################