文章标题
2016-11-06 22:39
281 查看
总体比例的区间估计
attach(faithful) #获取火山灰数据 population <- sample(eruptions,1000,replace = T) #获取总体 N <- length(population) #计算总体个数 p<-length(population[population>2.5])/N #计算population>2.5占总体的比例 Pi <- p #总体比例为Pi Pi ##[1] 0.658 layout(matrix(1:1, 1, 1)) n <- 30 z<-qnorm(0.975) #计算置信水平为95%的z值 plot(Pi,experimentes,type="n",xlim=c(Pi-6*sigma,Pi+6*sigma), ylim=c(0,experimentes)) abline(v=Pi) #做100次实验 experimentes=100 for(i in 1:experiments){ Pi_of_sample <- sample(population, n) #从总体随机取30个样本 Pi_of_x <- length(Pi_of_sample[Pi_of_sample>2.5])/n #计算大于2.5的比例 sigma = sqrt(Pi_of_x*(1-Pi_of_x)/n) co.inter <- c(Pi_of_x - z*sigma, Pi_of_x + z*sigma) #计算总体比例的置信区间 if(co.inter[1] < Pi & Pi <= co.inter[2]) lines(co.inter, c(i, i), type="l") else lines(co.inter, c(i, i), type="l", col=2) #如果总体比例不在该区间内, 用红线表示,否则,用黑线 } #做1000次实验 experiments=1000 success <- 0 for(i in 1:experiments){ Pi_of_sample <- sample(population, n) Pi_of_x <- length(Pi_of_sample[Pi_of_sample>2.5])/n sigma = sqrt(Pi_of_x*(1-Pi_of_x)/n) co.inter <- c(Pi_of_x - z*sigma, Pi_of_x + z*sigma) if(co.inter[1] < Pi & Pi <= co.inter[2]) success <- success +1 } success/experiments ![效果截图](http://img.blog.csdn.net/20161106223843396)
总体方差的区间估计
attach(faithful) population <- sample(eruptions,1000,replace = T) sigma <- var(population)*((N-1)/N) #总体方差1.306628 layout(matrix(1:1, 1, 1)) experimentes=100 n <- 30 kafang_left = qchisq(0.975,29) kafang_right = qchisq(0.025,29) plot(sigma,experimentes,type="n",xlim=c(sigma-3*sigma,sigma+3*sigma),ylim=c(0,experimentes)) abline(v=sigma) for(i in 1:experiments){ sigma_of_x = var(sample(population, n)) sigma_of_x co.inter <- c((n-1)*sigma_of_x/kafang_left, (n-1)*sigma_of_x/kafang_right) co.inter if(co.inter[1] < sigma & sigma <= co.inter[2]) lines(co.inter, c(i, i), type="l") else lines(co.inter, c(i, i), type="l", col=2) } experiments=100 success <- 0 for(i in 1:experiments){ sigma_of_x = var(sample(population, n)) co.inter <- c((n-1)*sigma_of_x/kafang_left, (n-1)*sigma_of_x/kafang_right) if(co.inter[1] < sigma & sigma <= co.inter[2]) success <- success +1 } success/experiments