rm(list=ls()) n <- 100 m <- 100 k <- 1000 alpha <- 0.05 p1 <- 0.5 p2 <- 0.3 hmus1 <- numeric(k) hmus2 <- numeric(k) for (i in 1:k){ x <- rbinom(n,1,prob=p1) #生成x服从二项分布随机数 hmus1[i]<- mean(x) } for (i in 1:k){ y <- rbinom(m,1,prob=p2) #生成y服从二项分布随机数 hmus2[i]<- mean(y) } hmu1 <- mean(hmus1) hmu2 <- mean(hmus2) hmu1 hmu2 sx2 <- sum((hmus1-hmu1)^2) sy2 <- sum((hmus2-hmu2)^2) sw2 <- (sx2+sy2)/(n+m+2) hup <- hmus1-hmus2+qt(alpha/2, n+m-2, lower.tail = FALSE)*sqrt(1/n+1/m)*sqrt(sw2) hlo <- hmus1-hmus2-qt(alpha/2, n+m-2, lower.tail = FALSE)*sqrt(1/n+1/m)*sqrt(sw2) hprob <- mean((hlo<=p1-p2)&(p1-p2<=hup)) hprob估计置信水平，哪里出错了

在计算 hprob 值时出现了错误。hprob 的计算方法应该是：计算差值（hlo<=p1-p2）与（p1-p2<=hup）同时成立的概率，即两者的交集，而不是计算两者的平均值。正确的计算方法是：

hprob <- mean(hlo<=p1-p2 & p1-p2<=hup)

这样计算得出的 hprob 才是估计置信水平的正确结果。

相关问题

置信区间报错rm(list=ls()) n <- 100 m <- 100 k <- 1000 alpha <- 0.05 p1 <- 0.5 p2 <- 0.3 hmus1 <- numeric(k) hmus2 <- numeric(k) for (i in 1:k){ x <- rbinom(n,1,prob=p1) #生成x服从二项分布随机数 hmus1[i]<- mean(x) } for (i in 1:k){ y <- rbinom(m,1,prob=p2) #生成y服从二项分布随机数 hmus2[i]<- mean(y) } hmu1 <- mean(hmus1) hmu2 <- mean(hmus2) hmu1 hmu2 sx2 <- sum((hmus1-hmu1)^2)/(k-1) sy2 <- sum((hmus2-hmu2)^2)/(k-1) sw2 <- ((n-1)*sx2+(m-1)*sy2)/(n+m-2) hup <- (hmus1-hmus2)+qt(1-alpha/2, df=n+m-2)*sqrt(1/n+1/m)*sqrt(sw2) hlo <- (hmus1-hmus2)-qt(1-alpha/2, df=n+m-2)*sqrt(1/n+1/m)*sqrt(sw2) hprob <- mean((hlo<=p1-p2) & (p1-p2<=hup)) hprob hp <- mean((hlo<=p1-p2)&(p1-p2<=hup)) hp

该段代码中的问题可能是在计算置信区间时使用了t分布的分位数，但是在计算方差时使用了样本方差而不是样本标准差的平方。应该修改为：

sx <- sqrt(sum((hmus1-hmu1)^2)/(k-1))
sy <- sqrt(sum((hmus2-hmu2)^2)/(k-1))
sw <- ((n-1)*sx^2+(m-1)*sy^2)/(n+m-2)
hup <- (hmus1-hmus2)+qt(1-alpha/2, df=n+m-2)*sqrt(1/n+1/m)*sqrt(sw)
hlo <- (hmus1-hmus2)-qt(1-alpha/2, df=n+m-2)*sqrt(1/n+1/m)*sqrt(sw)

这样就可以正确计算置信区间了。

rm(list=ls()) n <- 100 m <- 100 k <- 1000 alpha <- 0.05 p1 <- 0.5 p2 <- 0.3 hmus1 <- numeric(k) hmus2 <- numeric(k) for (i in 1:k){ x <- rbinom(n,1,prob=p1) #生成x服从二项分布随机数 hmus1[i]<- mean(x) } for (i in 1:k){ y <- rbinom(m,1,prob=p2) #生成y服从二项分布随机数 hmus2[i]<- mean(y) } hmu1 <- mean(hmus1) hmu2 <- mean(hmus2) hmu1 hmu2 sx2 <- sum((hmus1-hmu1)^2)/(k-1) sy2 <- sum((hmus2-hmu2)^2)/(k-1) sw2 <- ((n-1)*sx2+(m-1)*sy2)/(n+m-2) hup <- (hmus1-hmus2)+qt(1-alpha/2, df=n+m-2)*sqrt(1/n+1/m)*sqrt(sw2) hlo <- (hmus1-hmus2)-qt(1-alpha/2, df=n+m-2)*sqrt(1/n+1/m)*sqrt(sw2) hprob <- mean((hlo<=p1-p2) & (p1-p2<=hup)) hprob报错

这段代码的问题是在计算置信区间的时候，使用了错误的方差计算公式。应该修改为：

rm(list=ls()) 
n <- 100 
m <- 100 
k <- 1000 
alpha <- 0.05 
p1 <- 0.5 
p2 <- 0.3 
hmus1 <- numeric(k) 
hmus2 <- numeric(k) 

for (i in 1:k){   
  x <- rbinom(n,1,prob=p1)  #生成x服从二项分布随机数   
  hmus1[i]<- mean(x) 
}   

for (i in 1:k){   
  y <- rbinom(m,1,prob=p2)  #生成y服从二项分布随机数   
  hmus2[i]<- mean(y) 
}  

hmu1 <- mean(hmus1) 
hmu2 <- mean(hmus2) 
sx2 <- sum((hmus1-hmu1)^2)/(k-1) 
sy2 <- sum((hmus2-hmu2)^2)/(k-1) 
sw <- sqrt(((n-1)*sx2+(m-1)*sy2)/(n+m-2))  # 修改为sw
hup <- (hmus1-hmus2)+qt(1-alpha/2, df=n+m-2)*sqrt(1/n+1/m)*sw  
hlo <- (hmus1-hmus2)-qt(1-alpha/2, df=n+m-2)*sqrt(1/n+1/m)*sw   

hprob <- mean((hlo<=p1-p2) & (p1-p2<=hup)) 
hprob 

这样就可以正常运行并得到置信区间覆盖概率了。