R語言模擬置信區間估計
1.1總體均值的區間估計
方差已知,大樣本
程式碼:
attach(faithful) ##獲取火山灰資料
population <- sample(eruptions,1000,replace = T) #做1000次有放回取樣
N <- length(population) #總體數量
mu <- mean(population) #總體均值
#sd計算樣本方差(注:樣本方差除以(n-1))
sigma <- sd(population)*sqrt((N-1)/N) #總體標準差
layout(matrix(1:1, 1, 1))
experimentes=100
n <- 30
sample_sd <- sigma/sqrt(n)
z<-qnorm(0.975) #計算置信水平為95%的z值
z
plot(mu,experimentes,type="n",xlim=c(mu-6*sample_sd,mu+6*sample_sd),
ylim=c(0,experimentes))
abline(v=mu) #x=mu的垂直豎線作為參考線
##重複100次實驗,畫出每次的置信水平,若均值不在此區間,用紅線表示
for(i in 1:experiments){
mean_of_x <- mean(sample(population, n))
co.inter <- c(mean_of_x - z*sample_sd, mean_of_x + z*sample_sd)
if(co.inter[1] < mu & mu <= co.inter[2]) lines(co.inter, c(i, i), type="l")
else lines(co.inter, c(i, i), type="l", col=2)
}
#當樣本足夠大時,計算抽樣成功的次數所佔抽樣次數的比例,結果會接近95%
experiments=100000
success <- 0
for(i in 1:experiments){
mean_of_x <- mean(sample(population, n))
co.inter <- c(mean_of_x - z*sample_sd, mean_of_x + z*sample_sd)
if(co.inter[1] < mu & mu <= co.inter[2]) success <- success +1
}
success/experiments
效果截圖:
方差未知,小樣本
程式碼:
attach(faithful) ##獲取火山灰資料
population <- sample(eruptions,1000,replace = T) #做1000次有放回取樣
N <- length(population) #總體數量
mu <- mean(population) #總體均值
#sd計算樣本方差(注:樣本方差除以(n-1))
sigma <- sd(population)*sqrt((N-1)/N) #總體標準差
layout(matrix(1:1, 1, 1))
experimentes=100
n <- 9
tvalue<-qt(0.975,n-1) #計算置信水平為95%的t值
plot(mu,experimentes,type ="n",xlim=c(mu-6,mu+6),ylim=c(0,experimentes))
abline(v=mu)
##重複100次實驗,畫出每次的置信水平,若均值不在此區間,用紅線表示
for(i in 1:experiments){
samp <- sample(population, n)
mean_of_x <-mean(samp)
sd_of_x<-sd(samp)/sqrt(n)
co.inter <- c(mean_of_x - tvalue*sd_of_x, mean_of_x + tvalue*sd_of_x)
if(co.inter[1] < mu & mu <= co.inter[2]) lines(co.inter, c(i, i), type="l")
else lines(co.inter, c(i, i), type="l", col=2)
}
#當樣本足夠大時,計算抽樣成功的次數所佔抽樣次數的比例,結果會接近95%
experiments=100000
success <- 0
for(i in 1:experiments){
samp <- sample(population, n)
mean_of_x <-mean(samp)
sd_of_x<-sd(samp)/sqrt(n)
co.inter <- c(mean_of_x - tvalue*sd_of_x, mean_of_x + tvalue*sd_of_x)
if(co.inter[1] < mu & mu <= co.inter[2]) success <- success +1
}
success/experiments