(正規分布の)分散を推定する

  • 母分布が平均m、分散s^2の正規分布であるとき、そこからのランダムサンプルn個を取り、母分布のパラメタを推定する
m <- pi
s2 <- exp(1)
n <- 20
n.trial <- 10000
sample.mean <- sample.var <- sample.ub.var <- rep(0,n.trial)
for(i in 1:n.trial){
	X <- rnorm(n,m,sqrt(s2))
	sample.mean[i] <- mean(X)
	sample.var[i] <- sum((X-sample.mean[i])^2)/n
	sample.ub.var[i] <- var(X)
}
mean(sample.mean)
m
mean(sample.var)
mean(sample.ub.var)
s2
  • サンプリングを繰り返すと、標本平均のその平均は母分布の平均に近く、標本分散のその平均は母分布の分散からずれていて、標本から計算した不偏分散のその平均は母分布の分散に近いことが確かめられる
> mean(sample.mean)
[1] 3.148297
> m
[1] 3.141593
> mean(sample.var)
[1] 2.592012
> mean(sample.ub.var)
[1] 2.728434
> s2
[1] 2.718282
  • 全体として見れば、標本分散は標本からの不偏分散より常に小さく、その値は母分布の分散より小さいことも大きいこともある
matplot(cbind(sort(sample.var),sort(sample.ub.var)),type="l",cex=0.1)
abline(h=s2)

  • サンプリングを繰り返せればよいけれど、一度の標本だったとしたら、標本分散と標本からの不偏分散とでは、どちらが母分布の分散に近いのだろう
    • 差の平均は標本分散の方が小さく、差の分散も標本分散の方が小さい
diff.s.v <- abs(sample.var-s2)
diff.s.uv <- abs(sample.ub.var-s2)
mean(diff.s.v)
mean(diff.s.uv)
var(diff.s.v)
var(diff.s.uv)
> mean(diff.s.v)
[1] 0.6800217
> mean(diff.s.uv)
[1] 0.69536
> var(diff.s.v)
[1] 0.2519946
> var(diff.s.uv)
[1] 0.2905178
    • 差が標本分散の方が小さい場合と、逆に標本分散の方が大きい場合とに2分されている様子が次の図。y=xの赤線より下は標本分散の方が真からのずれが大きく、y=xの赤線より上は、標本分散の方が真に近い標本に対応する
plot(diff.s.v,diff.s.uv)
abline(0,1,col=2)
abline(h=0,col=3)
abline(v=0,col=3)

  • 標本分散と標本からの不偏分散とで母分散に近いのが標本分散である割合は以下のように0.5より少ないことから、標本からの不偏分散の方が母分布の分散をよりよく当てている確率は高いが、ときに、大きく外すこともある(差の平均が大きくなっている)ことがわかる
length(which(diff.s.v<diff.s.uv)) / length(diff.s.v)
> length(which(diff.s.v<diff.s.uv)) / length(diff.s.v)
[1] 0.4283
  • ずれの分布をみれば、差が大きい方場合が不偏分散の方に多いことがわかる
matplot(cbind(sort(diff.s.v),sort(diff.s.uv)),type="l",cex=0.1)

  • nの増加でどう変わるかをみる

m <- pi/4
s2 <- exp(1)
ns <- 2:30

mean.diffs <- var.diffs <- ratio.betters <- rep(0,length(ns))

for(nn in 1:length(ns)){
n <- ns[nn]
n.trial <- 10000
sample.mean <- sample.var <- sample.ub.var <- rep(0,n.trial)
for(i in 1:n.trial){
	#X <- sample(0:1,n,replace=TRUE,prob=c(1-m,m))
	X <- rnorm(n,m,sqrt(s2))
	sample.mean[i] <- mean(X)
	sample.var[i] <- sum((X-sample.mean[i])^2)/n
	sample.ub.var[i] <- var(X)
}
mean(sample.mean)
m
mean(sample.var)
mean(sample.ub.var)
s2
matplot(cbind(sort(sample.var),sort(sample.ub.var)),type="l",cex=0.1)
abline(h=s2)
diff.s.v <- abs(sample.var-s2)
diff.s.uv <- abs(sample.ub.var-s2)
mean(diff.s.v)
mean(diff.s.uv)
mean.diffs[nn] <- mean(diff.s.v)-mean(diff.s.uv)
var(diff.s.v)
var(diff.s.uv)
var.diffs[nn] <- var(diff.s.v)-var(diff.s.uv)

plot(diff.s.v,diff.s.uv)
abline(0,1,col=2)
abline(h=0,col=3)
abline(v=0,col=3)

ratio.betters[nn] <- length(which(diff.s.v<diff.s.uv)) / length(diff.s.v)
matplot(cbind(sort(diff.s.v),sort(diff.s.uv)),type="l",cex=0.1)
}

par(mfcol=c(2,2))
plot(ns,mean.diffs,type="b")
abline(h=0)
plot(ns,var.diffs,type="b")
abline(h=0)
plot(ns,ratio.betters,type="b")
abline(h=0.5)
par(mfcol=c(1,1))