2017-09-19 76 views
2

我從平均值爲10,標準差爲2的正態分佈中繪製了100個10號樣本。代碼如下:繪製從正態分佈得出的平均值的95%置信區間

n <- 10 
nreps<-100 
sample.mean<-numeric(nreps) 
for (i in 1:nreps) { 
    sample <- rnorm(n=n, mean = 10, sd = 2) 
    sample.mean[i] <- mean(sample) 
    a <- qnorm(0.95*2/sqrt(n)) 
    ci <- a 
} 
plot(sample.mean, 1:100) 

我希望創建一個看起來像這樣

這圖是我目前有

我知道我需要解釋每個平均值的左側和右側邊界,然後在它們之間插入一條水平線。超過95%置信區間的手段應該與其他手段的顏色不同。我剛剛開始學習R,所以非常感謝您的幫助。

回答

1

試試這樣說:

library(ggplot2) 
set.seed(1321) 
n <- 10 
sd <- 2 
n.reps <- 100 
my.mean <- 10 
alpha <- 0.05 

mydata <- matrix(rnorm(n = n.reps * n, mean = my.mean, sd =sd), ncol = n.reps) 

sample.means <- apply(mydata, 2, mean) 

error <- apply(mydata, 2, function(x) qt(p=1-alpha/2,df=length(x)-1)*sd (x)/sqrt(length(x))) 


dfx <- data.frame(sample.means, error, lcl = sample.means-error, ucl = sample.means+error, trial = 1:n.reps) 
dfx$miss <- dfx$ucl < my.mean | dfx$lcl > my.mean 
ggplot(dfx, aes(x = sample.means, y = trial, xmin = lcl, xmax = ucl, color = miss)) + geom_errorbarh() + geom_point(pch = 1) + 
    geom_vline(aes(xintercept=my.mean), lty=2) + 
    xlab("True Mean in Blue and 95% Confidence Intervals") + ylab ("Trial") + ggtitle(paste("Successful CI's:", 100*mean(!dfx$miss), "%")) + scale_color_manual(values = c("green", "red")) + 
    theme_bw() 

或使用基地:

oldpar <- par(xpd=FALSE) 
par(mar=c(8.1, 3.1, 3.1, 4.1)) 

with(subset(dfx, !miss), plot(sample.means, trial, 
           xlab = "Sample Mean", 
           ylab = "Trial", 
           col = "forestgreen", 
    xlim=c(min(dfx$lcl), max(dfx$ucl)))) 
with(subset(dfx, miss), points(sample.means, trial, 
           col = "red")) 

with(subset(dfx, miss), segments(lcl, trial, ucl, trial, col = "red")) 
with(subset(dfx, !miss), segments(lcl, trial, ucl, trial, col = "forestgreen")) 
abline(v = my.mean, lty = 2, lwd = 2, col = "blue") 

par(xpd=TRUE) 
legend("bottomright", c("Successful CI", "Miss"), lty = c(1,1), col = c("forestgreen", "red"), 
     inset=c(-0.1,-0.45)) 

title(main = paste("Successful CI's:", 100*mean(!dfx$miss), "%"), 
     sub = "True mean (in blue) and CI's") 
par(oldpar) 

HTH 詹姆斯

+0

如果你想使用qnorm正態分佈替代QT –

相關問題