標本平均の平均が母平均に一致することを表す例

「標本平均の平均が母平均に一致する」,これを説明するために,次のような例を示して教えてます。

  • 100個のデータを作って「これを母集団とする」とします。
  • その中からランダムに10個取り出して,標本平均を出します。
  • 標本を何度も取り出し,その度に標本平均を出します。
  • 標本平均の平均が母平均に近づいていくことを示します。

この説明をするために,100個の母集団を図示し,ランダムに抜き出す10個を色付けして表示,スライドに一枚ずつこの図を貼っていって,標本平均を計算,その平均を計算,というパターンを作ります。

これって誰得?と思うような図の書き方ですが,参考までに。

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()
library(ggrepel)
set.seed(20180531)
N <- 100
x <- round(rnorm(N, 50, 10))
matrix(x, ncol = 10)
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]   59   64   49   47   46   53   46   46   38    72
##  [2,]   49   51   63   58   52   59   58   40   64    45
##  [3,]   54   66   41   52   62   47   41   67   49    54
##  [4,]   51   67   66   53   62   41   32   76   60    54
##  [5,]   49   62   48   60   56   55   55   68   49    49
##  [6,]   42   51   65   47   44   30   54   48   60    41
##  [7,]   48   60   54   45   45   48   53   64   44    57
##  [8,]   55   61   44   49   66   33   62   65   32    49
##  [9,]   49   60   65   51   39   59   50   43   53    55
## [10,]   47   45   70   44   27   67   64   34   41    50
mean((x - mean(x))^2)/10
## [1] 9.83844
mean(x)
## [1] 52.34
dim(x) <- c(10, 10)
apply(x, 2, mean)
##  [1] 50.3 58.7 56.5 50.6 49.9 49.2 51.5 55.1 49.0 52.6
mean(apply(x,2,mean))
## [1] 52.34
d1 <- rep(1:10,10)
d2 <- rep(1:10,each=10)
dat <- transform(list(x=d1,y=d2,val=as.vector(x)))
ggplot(dat,aes(x=d1,y=d2,label=val))+geom_text(size=10)+xlab("")+ylab("")+
  theme(
    panel.background = element_blank(),
    panel.grid = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  )+scale_y_continuous(breaks=NULL)+scale_x_continuous(breaks=NULL)

# 何行何列目かを見つける関数
finds <- function(x,v){
  result <- c()
  for(i in 1:length(v)){
    pos <- which(x==v[i])
    if(length(pos)>1){pos=sample(pos,1)}
    result <- c(result,pos)
  }
  return(result)
}
 

# 可視化 ---------------------------------------------------------------------
n <- 10
X <- 10
samp.df <- data.frame(id=numeric(),val=numeric())
for(i in 1:X){
  dat$col <- 1
  s <- sample(x,n,replace = F)
  gr <- finds(x,s)
  dat[gr,]$col <- 2
  dat$col <- as.factor(dat$col)
  g <- ggplot(dat,aes(x=d1,y=d2,label=val,color=col))+geom_text(size=10)+xlab("")+ylab("")+
    theme(
      panel.background = element_blank(),
      panel.grid = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      legend.position = 'none'
    )+scale_y_continuous(breaks=NULL)+scale_x_continuous(breaks=NULL)+
    scale_color_manual(values = c("black","red"))
  plot(g)
  # 実際はプロットせずに図のファイルとして保存させます。s
  #ggsave(file = paste0("sample",i,".png"), plot = g, dpi = 100, width = 10,height = 5)
  samp.df <- rbind(samp.df,cbind(rep(i,n),s))
}

samp.df$V1 <- as.factor(samp.df$V1)
tapply(samp.df$s,samp.df$V1,mean)
##    1    2    3    4    5    6    7    8    9   10 
## 51.1 49.1 45.0 49.7 53.6 50.0 49.4 58.2 50.2 53.4
mean(samp.df$s)
## [1] 50.97
mean((tapply(samp.df$s,samp.df$V1,mean)-mean(x))^2)
## [1] 12.863
Share