不偏分散と標本分散の比較

不偏分散(\(N-1\)で割る方)が標本分散よりも分散の推定値として良いことを可視化します。

set.seed(20180704)
library(tidyverse)
## ─ Attaching packages ──────────────── tidyverse 1.2.1 ─
## ✔ ggplot2 3.0.0     ✔ purrr   0.2.5
## ✔ tibble  1.4.2     ✔ dplyr   0.7.6
## ✔ tidyr   0.8.1     ✔ stringr 1.3.1
## ✔ readr   1.1.1     ✔ forcats 0.3.0
## ─ Conflicts ───────────────── tidyverse_conflicts() ─
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
# 母集団のサイズ
N <- 100
# 乱数生成
x <- round(rnorm(N, 50, 10))
# 標本のサイズ
n <- 10
# 反復回数
X <- 100
samp.df <- data.frame(id=numeric(),val=numeric())
# サンプリングを繰り返す
for(i in 1:X){
  s <- sample(x,n,replace = F)
  samp.df <- rbind(samp.df,cbind(rep(i,n),s))
}
samp.df$V1 <- as.factor(samp.df$V1)

# 標本分散
bV <- function(x){
  mx <- mean(x)
  v <- mean((x-mx)^2)
  return(v)
}

# サンプリングごとの標本分散
mv <- tapply(samp.df$s,samp.df$V1,bV)

# 標本分散の平均ベクトル
mm <- c()
for(i in 1:length(mv)){
  tmp <- mean(mv[1:i])
  mm <- c(mm,tmp)
}

# サンプリングごとの不偏分散
mv2 <- tapply(samp.df$s,samp.df$V1,var)

# 不偏分散の平均ベクトル
mm2 <- c()
for(i in 1:length(mv2)){
  tmp <- mean(mv2[1:i])
  mm2 <- c(mm2,tmp)
}


# 可視化
cbind(mm,mm2) %>% transform() %>% mutate(rep=as.numeric(rownames(.))) %>% tidyr::gather(key,val,-rep,factor_key=T) %>%
  ggplot(aes(x=rep,y=val,color=key))+ geom_point()+geom_line()+geom_hline(yintercept=bV(x),color="red")+
  ylab('means of sample variance')+xlab('number of samples')+theme_grey(base_family = "HiraKakuProN-W3")+
  scale_colour_discrete(name='分散の種類',labels=c("標本分散","不偏分散"))

Share