2017-09-19 19 views
2

私は平均10と標準偏差2の正規分布からサイズ10の100サンプルを抽出しました。以下のコード:平均分布から導き出された平均の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を使用したい場合 –

関連する問題