標本抽出法

立命館大学産業社会学部教授 筒井淳也 様 のスライド で学ばせていただいた記録です。標本抽出の理論とシミュレーションについて理解を深めることができました。御礼申し上げます。

1 母集団の準備.

1クラス40人、全体は10クラスで400人のある学年。 番号の若いクラスほど、成績が悪い傾向がある。
40人を標本調査することで、この学年のテストの点の平均（＝母数）を知りたい。
母集団を固定する.

set.seed(1)#乱数の種を指定して母集団を固定する
d<-list (10) #リストのi番目は40個の乱数
for (i in 1:10){
  d[[i]]<-rnorm(n=40, mean=50,sd=10)+3*i-16.5
}
  d2<-c(d[[1]],d[[2]],d[[3]],d[[4]],d[[5]],d[[6]],d[[7]],d[[8]],d[[9]],d[[10]])
head(d2)#40個の乱数10セットで合計400個

## [1] 30.23546 38.33643 28.14371 52.45281 39.79508 28.29532

par(family = "HiraKakuProN-W3") 
boxplot (d[[1]],d[[2]],d[[3]],d[[4]],d[[5]],d[[6]]
                  ,d[[7]],d[[8]],d[[9]],d[[10]])

#axis (1,seq(1:100))
mean(d2);var(d2)*399/400#推定したい母平均

## [1] 50.38089

## [1] 166.813

1.0.1 群間変動.

SSB=c()
for (i in 1:10){
  ssb <- 40*(mean(d[[i]])-mean(d2))^2
  SSB <-c(SSB,ssb)
}
#群間分散
sum(SSB)/400

## [1] 74.88577

1.0.2 群内変動.

SSW=c()
for (i in 1:10){
  ssw <- sum((d[[i]]-mean(d[[i]]))^2)
  SSW <-c(SSW,ssw)
}
#群内分散
sum(SSW)/400

## [1] 91.92727

1.0.3 全変動.

SST <- sum(SSB)+sum(SSW)
#母分散
SST/400

## [1] 166.813

#40個の標本を抽出したときの標本分散
SST/400/40

## [1] 4.170326

2 単純無作為抽出.

 m1<-numeric (10000)#10000個のデータの枠を作る
#400個から40個を標本抽出しその平均を出す。この操作を10000回行う
for(i in 1:10000){
  m1[i]<-mean (sample(x=d2,size=40, replace=T))
}
head(m1)

## [1] 53.78538 50.65396 49.13630 48.45420 49.54169 51.61348

2.1 推定値の分布.

400人から40人を抽出した標本平均の分布.

par(family = "HiraKakuProN-W3") 
hist(m1,xlim=c(35,65),xlab="平均",ylab="回数",main="単純無作為抽出")
#母数
abline (v=mean (d2))

mean (m1)

## [1] 50.35486

var (m1)

## [1] 4.193911

var(d2)*399/400/40#標本の平均の分散の理論値

## [1] 4.170326

3 多段抽出.

m2<-numeric (10000)#10000個のデータ
#1から10の一様乱数を用いて d[[]]で4つのリストを選び、各リストの40個から10個選ぶ
for (i in 1:10000) {
m2[i]<-mean(c(
  sample(x=d[[floor(runif(n=1,min=1,max=11))]],replace=T,size=10),
  sample(x=d[[floor(runif(n=1,min=1,max=11))]],replace=T,size=10),
  sample(x=d[[floor(runif(n=1, min=1, max=11))]],replace=T,size=10),
  sample(x=d[[floor(runif(n=1, min=1, max=11))]],replace=T,size=10)
))
}
str(m2)#多段平均40個を10000セット実施しヒストグラム

##  num [1:10000] 55.1 60.3 51 46.9 41.9 ...

par(family = "HiraKakuProN-W3") 
hist(m2,xlim=c(35,65),xlab="平均",ylab="回数",main="多段抽出")
abline (v=mean(d2))

mean (m2)

## [1] 50.31489

var (m2)

## [1] 21.32038

3.1 理論値.

Vts <-sum(SSW)/400/40+sum(SSB)/400/4
Vts

## [1] 21.01963

多段抽出では分散が大きくなる.

4 層化抽出法.

各層からまんべんなく４個ずつ計40個を抽出する。

m3<-numeric(10000)
for (i in 1:10000) {
  m3[i]<-mean (c(
             sample(x=d[[1]], replace=T,size=4),
             sample(x=d[[2]], replace=T,size=4),                 
             sample(x=d[[3]], replace=T,size=4),
             sample(x=d[[4]], replace=T,size=4),
             sample(x=d[[5]], replace=T,size=4),
             sample(x=d[[6]], replace=T,size=4),
             sample(x=d[[7]], replace=T,size=4),
             sample(x=d[[8]], replace=T,size=4),
             sample(x=d[[9]], replace=T,size=4),
             sample(x=d[[10]], replace=T,size=4)
             ))
}
var (m3)

## [1] 2.280152

4.1 理論値.

sum(SSW)/400/40

## [1] 2.298182

群内変動から求めた分散に一致する。

5 比較.

母集団の平均と分散が \((\mu,\sigma^2)\) であるとする。
母集団は M 個の群からなり、群 s の平均と分散は \((\mu_s,\sigma_s^2)\) とする。
　　　\(\displaystyle \sigma^2=\sum_{s=1}^M W_s\sigma_s^2+ \sum_{s=1}^MW_s(\mu_s-\mu)^2\)　（郡内分散＋群間分散）　

5.1 標本分散.

・層別抽出では、群内分散を標本数で割ったものになる。（小）.
　　　\(\displaystyle \dfrac{1}{n}\sum_{s=1}^M W_s\sigma_s^2\).
・無作為抽出では、郡内分散と群間分散を標本数で割ったものになる。（中）.
　　　\(\displaystyle \dfrac{1}{n}\sum_{s=1}^M W_s\sigma_s^2+\dfrac{1}{n} \sum_{s=1}^MW_s(\mu_s-\mu)^2\)　

・２段抽出では、群間分散は標本数で割らない（大）.
　　　\(\displaystyle \dfrac{1}{n}\sum_{s=1}^M W_s\sigma_s^2+ \sum_{s=1}^MW_s(\mu_s-\mu)^2\)　

5.2 具体例.

par (mfrow=c (3,1), oma=c(0,0,0,0))#omaは余白調整
par(family = "HiraKakuProN-W3") 
hist(m1,xlim=c(35,65),xlab="平均",ylab="回数",main="単純無作為抽出")
abline(v=mean (d2))

hist(m2,xlim=c(35,65),xlab="平均",ylab="回数",main="多段抽出")
abline(v=mean (d2))

hist(m3,xlim=c(35,65),xlab="平均",ylab="回数",main="層化抽出")
abline(v=mean (d2))

多段抽出＞単純無作為抽出＞層化抽出.

5.3 層化抽出の優位性.

m4 <- numeric (10000)
#復元抽出
for(i in 1:10000){
m4[i]<-mean (sample (x=d2 ,size=40*1.8, replace=T)) 
}
var(m1)#単純無作為抽出

## [1] 4.193911

var(m4)#標本1.8倍の単純無作為抽出

## [1] 2.340814

var(m3)#層化抽出

## [1] 2.280152

分散を1/1.8にするには標本サイズnが1.8倍必要.
層化抽出が標本サイズを約半分に節約している

5.4 層化抽出が有効である条件.

①母集団における各層(グループ）の割合がわかっている.
　例:ある学年における各クラスの人数.
②標本抽出枠(名簿)に層のデータが含まれる.
　例：学年名簿に所属クラスが記されている.
③層ごとに母数が異なる.
　例:クラスごとにテストの平均得点が異なる.

①または②が満たされないとき、次善の策として単純無作為抽出がとられる。
③が満たされないとき、層化抽出の精度は単純無作為抽出と同じ.

5.5 標本抽出法を比べる.

	推定値の平均	推定値の分散	調査コスト
单純無作為抽出	○母数と同じ	△計算しやすい	△調査対象が散らばる
多段抽出	△割当次第で○	×層ごとに異なるほど大きくなる	○
層化抽出	△割当次第で○	○層ごとに異なるほど小さくなる	×調査対象が散らばる

5.6 層化抽出が優位である理由.

単純無作為抽出の場合
　① 標本平均の分散=母分散／標本サイズ.
　② 母分散=母集団の偏差平方和／母集団のサイズ.
　③ 母集団の偏差平方和=郡内平方和＋群間平方和（＃平方和分解）.

①に②③を代入すると
　④ 標本平均の分散＝(郡内平方和＋群間平方和）/母集団のサイズ/標本サイズ.

層化抽出の場合.
　⑤ 標本平均の分散＝Σ(i組の分散✕i組のサイズ)／母集団のサイズ／標本サイズ.
　⑥ 郡内平方和=Σ（i組の分散×i組のサイズ）.
⑤に⑥を代入すると.
　⑦ 標本平均の分散=群内平方和／母集団のサイズ／標本サイズ.

単純無作為抽出と層化抽出の差.
標本平均の分散の差=④ー⑦=群間平方和／母集団のサイズ／標本サイズ.

層ごとに平均（母数)が異なるほど、層化抽出の効果が大きい.

5.7 層化抽出優位の具体例.

#平方和を求める関数
sv <- function(x){
  return(var(x)*(length(x)-1))
}

#群間平方和
#定義した関数を利用
#群間変動
gk <- sv(c(
mean(d[[1]]),mean(d[[2]]),mean(d[[3]]),mean(d[[4]]),mean(d[[5]]),
mean(d[[6]]),mean(d[[7]]),mean(d[[8]]),mean(d[[9]]),mean(d[[10]])))*40
gk

## [1] 29954.31

#郡内平方和
gn <- sum(c(
sv(d[[1]]),sv(d[[2]]),sv(d[[3]]),sv(d[[4]]),sv(d[[5]]),
sv(d[[6]]),sv(d[[7]]),sv(d[[8]]),sv(d[[9]]),sv(d[[10]])))
gn

## [1] 36770.91

sv(d2)

## [1] 66725.22

SST=SSB+SSW が成立している。
400で割れば母分散 \(\sigma^2\) 、さらに40で割れば標本の分散となる。
次のvar(m1)で確かめられる。

#単純無作為抽出の標本平均の分散
var(m1)

## [1] 4.193911

#（郡内平方和＋群間平方和）/母集団のサイズ/標本サイズ
(gn+gk)/400/40

## [1] 4.170326

#層化抽出の標本平均の分散
var(m3)

## [1] 2.280152

#郡内平方和/母集団のサイズ/標本サイズ
gn/400/40

## [1] 2.298182

#単純無作為抽出と層化抽出の分散の差
var(m1)-var(m3)

## [1] 1.913759

#群内平方和/母集団のサイズ/標
gk/400/40

## [1] 1.872144

5.8 まとめ.

1.単純無作為抽出が常に最善ではない.
　既知の情報を意図的に揃えれば、小さいデータで精度のよい推定が可能.
2.普通の検定、推定は単純無作為抽出された標本の話.
　多段抽出、層化抽出された標本の場合の計算は複雑（誤差の独立性）.
3.Rで簡単にシミュレーションができる 　数式にもとづく理論値の正しさを確認できる.

6 配分法.

統計検定準１級ワーク_問21.1を題材にシミュレーションを行う。
## データリストの作成.

d<-list (4) 
d

## [[1]]
## [1] 4

d[[1]]<-rnorm(n=20, mean=15,sd=20)
d[[2]]<-rnorm(n=10, mean=150,sd=70)
d[[3]]<-rnorm(n=5, mean=510,sd=290)
d[[4]]<-rnorm(n=5, mean=1010,sd=80)
d2<-c(d[[1]],d[[2]],d[[3]],d[[4]])
d2

##  [1]   -5.8509149   14.9266012    3.3056357   24.7968767    3.8225032
##  [6]   41.0766017   43.2899601   33.1818854  -25.4334670   -1.4366323
## [11]   39.1203110   45.6461799    0.2099484   11.3485793    0.2696753
## [16]   23.1368022   14.2698601  -10.3527062   15.5795727   55.5218421
## [21]   50.2175633  210.5207058  108.5014259   99.5545915   62.7327863
## [26]  275.2851497  174.2007511  185.9845463  204.3007697  101.2985308
## [31] 1308.7925378  555.0826296  685.7876987   69.9029139  -75.3493853
## [36] 1177.2817405 1073.3918776  943.9514943  916.6895691 1037.9192864

sum(d2)

## [1] 9492.476

boxplot (d[[1]],d[[2]],d[[3]],d[[4]])

6.1 層化抽出.

各層から非復元抽出を行う.
各群からのサンプル数の決め方は、

• 等配分法
• 比例配分法（各層の大きさに比例）
• ネイマン配分法（各層の大きさと分散の大きさに比例）

がある。

6.1.1 等配分法.

n1=2
n2=2
n3=2
n4=2
m1<-numeric(10000)
for (i in 1:10000) {
  S1 <- sample(x=d[[1]], replace=F,size=n1)
  S2 <- sample(x=d[[2]], replace=F,size=n2)               
  S3 <- sample(x=d[[3]], replace=F,size=n3)
  S4 <- sample(x=d[[4]], replace=F,size=n4)
  Yhat <- 20*mean(S1)+10*mean(S2)+5*mean(S3)+5*mean(S4)
  m1[i]<-Yhat
}
hist(m1/40)

mean(m1/40)

## [1] 237.2678

var(m1/40)

## [1] 1663.511

sd(m1/40)

## [1] 40.78616

等配分法による推定値は、\(237\pm 40\) となる。

6.1.2 比例配分法

n1=4
n2=2
n3=1
n4=1
m2<-numeric(10000)
for (i in 1:10000) {
  S1 <- sample(x=d[[1]], replace=F,size=n1)
  S2 <- sample(x=d[[2]], replace=F,size=n2)               
  S3 <- sample(x=d[[3]], replace=F,size=n3)
  S4 <- sample(x=d[[4]], replace=F,size=n4)
  Yhat <- 20*mean(S1)+10*mean(S2)+5*mean(S3)+5*mean(S4)
  m2[i]<-Yhat
}
hist(m2/40)

mean(m2/40)

## [1] 237.0307

var(m2/40)

## [1] 4094

sd(m2/40)

## [1] 63.98437

比例配分法による推定値は、\(237\pm 64\) となる。
大きな分散の群の影響を受け、あまり良い分布とは言えない。

n1=1
n2=2
n3=4
n4=1
m3<-numeric(10000)
for (i in 1:10000) {
  S1 <- sample(x=d[[1]], replace=F,size=n1)
  S2 <- sample(x=d[[2]], replace=F,size=n2)               
  S3 <- sample(x=d[[3]], replace=F,size=n3)
  S4 <- sample(x=d[[4]], replace=F,size=n4)
  Yhat <- 20*mean(S1)+10*mean(S2)+5*mean(S3)+5*mean(S4)
  m3[i]<-Yhat
}
hist(m3/40)

mean(m3/40)

## [1] 237.4658

var(m3/40)

## [1] 613.7483

sd(m3/40)

## [1] 24.77394

比例配分法による推定値は、\(237\pm 25\) となる。
分布の形状も良く、精度の良い推定となっている。

復元抽出にすると次のようになる。

n1=1
n2=2
n3=4
n4=1
m3<-numeric(10000)
for (i in 1:10000) {
  S1 <- sample(x=d[[1]], replace=T,size=n1)
  S2 <- sample(x=d[[2]], replace=T,size=n2)               
  S3 <- sample(x=d[[3]], replace=T,size=n3)
  S4 <- sample(x=d[[4]], replace=T,size=n4)
  Yhat <- 20*mean(S1)+10*mean(S2)+5*mean(S3)+5*mean(S4)
  m3[i]<-Yhat
}
hist(m3/40)

mean(m3/40)

## [1] 236.8604

var(m3/40)

## [1] 1331.308

sd(m3/40)

## [1] 36.4871

精度は落ちるようである。