スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

標本分散と不偏分散

標本数$N$の標本分散が不偏分散の$\frac{N-1}{N}$倍となることの、elephantな「力技」による確認を行う。


標本分散は次の関数として定義される:
variance <- function(x){
mean(x^2 - mean(x)^2)
}

Rでvar()が返す値は不偏分散であり、var()=$\frac{N}{N-1}$variance()であることに注意する。


正規分布の場合:
pdf("gauss17Mar2013.pdf", height=5.5, width=7, pointsize=14)
x <- rnorm(100000,0,1)
hist(x,breaks=100,probability=T,main="Norm",col="#CDBEE9")
dev.off()
windows()

gauss17Mar2013_01-mini.png

test <- matrix(nrow=100000, ncol=200)
test2 <- matrix(nrow=100000, ncol=200)
plot <- matrix(nrow=200,ncol=1)
plot2 <- matrix(nrow=200,ncol=1)

for(j in c(5:200)){
for(i in c(1:100000)){
x <- rnorm(j,0,2)
test[i,j] <- variance(x)
test2[i,j] <- var(x)
}
plot[j,1] <- mean(test[,j])
plot2[j,1] <- mean(test2[,j])
}

pdf("normplot17Mar2013.pdf", height=5.5, width=7, pointsize=14)
plot(plot2[,1],xlim=c(4,200),ylim=c(3.1,4.1),xlab="N",ylab="variance",
main="Norm",col="#E82202",pch=20,ps=0.01)
par(new=T)
plot(plot[,1],xlim=c(4,200),ylim=c(3.1,4.1),xlab="N",ylab="variance",
col="#6013FB",pch=20,ps=0.01)
dev.off()
windows()

normplot17Mar2013_01-mini.png



Weibull分布の場合(※1)
pdf("weibull17Mar2013.pdf", height=5.5, width=7, pointsize=14)
x <- rweibull(100000,1,7)
hist(x,breaks=100,xlim=c(0,60),probability=T,main="Weibull",col="#CDBEE9")
dev.off()
windows()

weibull17Mar2013_01-mini.png

※1 なぜここでWeibull分布? と思われた方も多いだろう。特に理由は無いのだが、私はWeibull分布とtangledな関係にあるらしく、例に示すなら正規分布の次にWeibull分布、という約束に、今日からした。

for(j in c(5:200)){
for(i in c(1:100000)){
x <- rnorm(j,0,2)
test[i,j] <- variance(x)
test2[i,j] <- var(x)
}
plot[j,1] <- mean(test[,j])
plot2[j,1] <- mean(test2[,j])
}

pdf("normplot17Mar2013.pdf", height=5.5, width=7, pointsize=14)
plot(plot2[,1],xlim=c(4,200),ylim=c(3.1,4.1),xlab="N",ylab="variance",
main="Norm",col="#E82202",pch=20,ps=0.01)
par(new=T)
plot(plot[,1],xlim=c(4,200),ylim=c(3.1,4.1),xlab="N",ylab="variance",
col="#6013FB",pch=20,ps=0.01)
dev.off()
windows()

weibullplot17Mar2013_02-mini.png
スポンサーサイト
Profile

zoa

Author:zoa
Metropolis (1927)

Calendar
02 | 2013/03 | 04
- - - - - 1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30
31 - - - - - -
Labels
Admin
Previous Posts
Recent Comments
Recent Trackback
Archive
このページのトップへ
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。