用R语言求置信区间


用R语言求置信区间是很方便的，而且很灵活，至少我觉得比spss好多了。


如果你要求的只是95%的置信度的话，那么用一个很简单的命令就可以实现了

首先，输入da=c(你的数据，用英文逗号分割），然后t.test(da)，运行就能得到结果了。


我的数据是newbomb <-
c(28,26,33,24,34,-44,27,16,40,-2,29,22,24,21,25,30,23,29,31,19)


t.test(newbomb)得到的结果如下

 

如果要求任意置信度下的置信区间的话，就需要自己编一个函数了。

当然，有两点要记住的，置信区间的计算在知道方差和不知道方差的情况下，计算公式是不一样的。

下面做一个两种情况下都可以用的函数。

confint<-function(x,sigma=-1,alpha=0.05)

  {

     
n<-length(x)

     
xb<-mean(x)

     
if(sigma>=0)

     
    {

     
     
 tmp<-sigma/sqrt(n)*qnorm(1-alpha/2);df<-n


     
     }

     
 else{

     
   
 tmp<-sd(x)/sqrt(n)*qt(1-alpha/2,n-1);df<-
n-1

     
     }

     
 data.frame(mean=xb,df=df,a=xb-tmp,b=xb+tmp)

   }

这个函数的使用：

如果不知道方差，则confint（x,alpha）
 知道方差，则confint（x,sigma,alpha)

这样就能计算出结果了。http://www.cda.cn/view/946.html
							
		
转载于:https://www.cnblogs.com/amengduo/p/9587628.html

       最近工作中用到了区间估计，将自己程序中用到的方法整理如下，仅供参考~

interval_estimated <- function(x, sigma=-1, side=0, alpha=0.05){
  n <- length(x); xb <- mean(x)
  if (sigma >= 0){
    if (side < 0){
	  tmp <- sigma/sqrt(n)*qnorm(1-alpha)
	  a <- -Inf; b <- xb+tmp
	} else if (side > 0){
	  tmp <- sigma/sqrt(n)*qnorm(1-alpha)
	  a <- xb-tmp; b <- Inf
	} else{
      tmp <- sigma/sqrt(n)*qnorm(1-alpha/2)
	  a <- xb-tmp; b <- xb+tmp
    }
	df <- n
  } else{
    if (side < 0){
	  tmp <- sd(x)/sqrt(n)*qt(1-alpha,n-1)
	  a <- -Inf; b <- xb+tmp
	} else if(side > 0){
	  tmp <- sd(x)/sqrt(n)*qt(1-alpha,n-1)
	  a <- xb-tmp; b <- Inf
	} else{
	  tmp <- sd(x)/sqrt(n)*qt(1-alpha/2,n-1)
	  a <- xb-tmp; b <- xb+tmp
	}
	df <- n-1
  }
  data.frame(a=a, b=b)
}

       在上述方法中，x是数据构成的向量。sigma是总体的标准差，sigma的默认值为-1，当标准差已知的时候，输入相应的值，程序采用正态分布估计区间端点，否则，程序采用t分布估计区间端点。side参数控制单双侧置信区间，默认值为0，当side<0时，求单侧置信区间的上限，当side>0时，求单侧置信区间下限，当side=0时，求双侧置信区间。该方法输出为区间的两个端点的值。

$\mu = 5,\theta = 4,n = 10, K = 100, \alpha = 0.1$
CI_L<-NULL 
CI_U<-NULL
for (k in 1:100){
  samps<-rnorm(10,mean=5,sd=2) 
  mu<-mean(samps) 
  sd<-sd(samps) 
  critv<-qt(1-0.1/2,9)
  CI_L<-c(CI_L,mu-critv*sd/sqrt(10)) 
  CI_U<-c(CI_U,mu+critv*sd/sqrt(10))
}
plot(0:100,
     type='n',
     ylim =c(round(min(CI_L))-1,round(max(CI_U))+1), 
     xlab ='',
     ylab ='90% CI')
segments(1:100,CI_L,1:100,CI_U,col='red',lwd=2) 
abline(h=5)


实验内容一
1.
1.固定样本量$n = 100$和$\alpha = 0.05$,观察重复次数100、200和400时置信区间包含真值$\mu = 15$的频率是否接近置信度$1- \alpha = 0.95$
# k = 100
CI_L<-NULL 
CI_U<-NULL
for (k in 1:100){
  samps<-rnorm(100,mean=15,sd=2) 
  mu<-mean(samps) 
  sd<-sd(samps) 
  critv<-qt(1-0.05/2,99)
  CI_L<-c(CI_L,mu-critv*sd/sqrt(100)) 
  CI_U<-c(CI_U,mu+critv*sd/sqrt(100))
}
plot(0:100,
     type='n',
     ylim =c(round(min(CI_L))-1,round(max(CI_U))+1), 
     xlab ='',
     ylab ='95% CI')
segments(1:100,CI_L,1:100,CI_U,col='red',lwd=2) 
abline(h=15)


# k = 200
CI_L<-NULL 
CI_U<-NULL
for (k in 1:200){
  samps<-rnorm(100,mean=15,sd=2) 
  mu<-mean(samps) 
  sd<-sd(samps) 
  critv<-qt(1-0.05/2,99)
  CI_L<-c(CI_L,mu-critv*sd/sqrt(100)) 
  CI_U<-c(CI_U,mu+critv*sd/sqrt(100))
}
plot(0:200,
     type='n',
     ylim =c(round(min(CI_L))-1,round(max(CI_U))+1), 
     xlab ='',
     ylab ='95% CI')
segments(1:200,CI_L,1:200,CI_U,col='red',lwd=2) 
abline(h=15)


# k = 400
CI_L<-NULL 
CI_U<-NULL
for (k in 1:400){
  samps<-rnorm(100,mean=15,sd=2) 
  mu<-mean(samps) 
  sd<-sd(samps) 
  critv<-qt(1-0.05/2,99)
  CI_L<-c(CI_L,mu-critv*sd/sqrt(100)) 
  CI_U<-c(CI_U,mu+critv*sd/sqrt(100))
}
plot(0:400,
     type='n',
     ylim =c(round(min(CI_L))-1,round(max(CI_U))+1), 
     xlab ='',
     ylab ='95% CI')
segments(1:400,CI_L,1:400,CI_U,col='red',lwd=2) 
abline(h=15)


2.
2.设置$\alpha = 0.10$，其他保持1不变，重复1，观察模拟结果；并观察与1中置信区间长度对比效果（随的变化）
# k = 100
CI_L<-NULL 
CI_U<-NULL
for (k in 1:100){
  samps<-rnorm(10,mean=15,sd=2) 
  mu<-mean(samps) 
  sd<-sd(samps) 
  critv<-qt(1-0.1/2,90)
  CI_L<-c(CI_L,mu-critv*sd/sqrt(10)) 
  CI_U<-c(CI_U,mu+critv*sd/sqrt(10))
}
plot(0:100,
     type='n',
     ylim =c(round(min(CI_L))-1,round(max(CI_U))+1), 
     xlab ='',
     ylab ='90% CI')
segments(1:100,CI_L,1:100,CI_U,col='red',lwd=2) 
abline(h=15)


# k = 200
CI_L<-NULL 
CI_U<-NULL
for (k in 1:200){
  samps<-rnorm(100,mean=15,sd=2) 
  mu<-mean(samps) 
  sd<-sd(samps) 
  critv<-qt(1-0.1/2,99)
  CI_L<-c(CI_L,mu-critv*sd/sqrt(100)) 
  CI_U<-c(CI_U,mu+critv*sd/sqrt(100))
}
plot(0:200,
     type='n',
     ylim =c(round(min(CI_L))-1,round(max(CI_U))+1), 
     xlab ='',
     ylab ='90% CI')
segments(1:200,CI_L,1:200,CI_U,col='red',lwd=2) 
abline(h=15)


# k = 400
CI_L<-NULL 
CI_U<-NULL
for (k in 1:400){
  samps<-rnorm(10,mean=15,sd=2) 
  mu<-mean(samps) 
  sd<-sd(samps) 
  critv<-qt(1-0.1/2,90)
  CI_L<-c(CI_L,mu-critv*sd/sqrt(10)) 
  CI_U<-c(CI_U,mu+critv*sd/sqrt(10))
}
plot(0:400,
     type='n',
     ylim =c(round(min(CI_L))-1,round(max(CI_U))+1), 
     xlab ='',
     ylab ='90% CI')
segments(1:400,CI_L,1:400,CI_U,col='red',lwd=2) 
abline(h=15)


3.
3.将1中样本量变成n = 200,其他不变，重复1，观察模拟结果，并观察与1中置信区间长度对比效果（随n的变化）
# k = 100
CI_L<-NULL 
CI_U<-NULL
for (k in 1:100){
  samps<-rnorm(200,mean=15,sd=2) 
  mu<-mean(samps) 
  sd<-sd(samps) 
  critv<-qt(1-0.05/2,199)
  CI_L<-c(CI_L,mu-critv*sd/sqrt(200)) 
  CI_U<-c(CI_U,mu+critv*sd/sqrt(200))
}
plot(0:100,
     type='n',
     ylim =c(round(min(CI_L))-1,round(max(CI_U))+1), 
     xlab ='',
     ylab ='95% CI')
segments(1:100,CI_L,1:100,CI_U,col='red',lwd=2) 
abline(h=15)


# k = 200
CI_L<-NULL 
CI_U<-NULL
for (k in 1:200){
  samps<-rnorm(200,mean=15,sd=2) 
  mu<-mean(samps) 
  sd<-sd(samps) 
  critv<-qt(1-0.05/2,199)
  CI_L<-c(CI_L,mu-critv*sd/sqrt(200)) 
  CI_U<-c(CI_U,mu+critv*sd/sqrt(200))
}
plot(0:200,
     type='n',
     ylim =c(round(min(CI_L))-1,round(max(CI_U))+1), 
     xlab ='',
     ylab ='95% CI')
segments(1:200,CI_L,1:200,CI_U,col='red',lwd=2) 
abline(h=15)


# k = 400
CI_L<-NULL 
CI_U<-NULL
for (k in 1:400){
  samps<-rnorm(200,mean=15,sd=2) 
  mu<-mean(samps) 
  sd<-sd(samps) 
  critv<-qt(1-0.05/2,199)
  CI_L<-c(CI_L,mu-critv*sd/sqrt(200)) 
  CI_U<-c(CI_U,mu+critv*sd/sqrt(200))
}
plot(0:400,
     type='n',
     ylim =c(round(min(CI_L))-1,round(max(CI_U))+1), 
     xlab ='',
     ylab ='95% CI')
segments(1:400,CI_L,1:400,CI_U,col='red',lwd=2) 
abline(h=15)

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