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(1-alpha/2, n+m-2, lower.tail = FALSE)*sqrt(1/n+1/m)*sqrt(sw2) hlo <- hmus1-hmus2-qt(1-alpha/2, n+m-2, lower.tail = FALSE)*sqrt(1/n+1/m)*sqrt(sw2) hprob <- mean((hlo<=p1-p2) & (p1-p2<=hup)) hprob报错

在计算 hup 和 hlo 值时，应该使用 qt 函数的第二个参数 df 为自由度，而不是 n+m-2。正确的计算方法是：

hup <- hmus1 - hmus2 + qt(1-alpha/2, df=n+m-2, lower.tail=FALSE)*sqrt(1/n+1/m)*sqrt(sw2)
hlo <- hmus1 - hmus2 - qt(1-alpha/2, df=n+m-2, lower.tail=FALSE)*sqrt(1/n+1/m)*sqrt(sw2)

此外，需要注意代码中变量名的大小写是否一致，比如 hprob 和 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 <- ((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)*sqrt(sw) 
hlo <- (hmus1-hmus2)-qt(1-alpha/2, df=n+m-2)*sqrt(1/n+1/m)*sqrt(sw)  

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

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