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报错
时间: 2024-03-05 08:53:49 浏览: 17
这段代码的问题是在计算置信区间的时候,使用了错误的方差计算公式。应该修改为:
```
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
```
这样就可以正常运行并得到置信区间覆盖概率了。