1. 程式人生 > >R語言模擬置信區間估計

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