『R』により標本平均の信頼区間を求めて図示する方法

以下は、『R』により、標準正規分布（平均$\mu=0$、分散$\sigma^2=1$）から採取した200個の標本の平均値の信頼区間を100回求め、図示するスクリプトである。Wikipediaの「信頼区間」の項（リンク）に掲載されたような図を描画してみたかったので、信頼区間の計算式については、イチから作成した。ただし、基本パッケージに収録された主要な統計的検定手法では、信頼区間も併せて出力されるし、bootパッケージやnlmeパッケージには、それぞれのパッケージに則した（ブートストラップ法で信頼区間を求める方法やモデルのパラメータについての信頼区間を表示する）方法が収録されている。このため、フルスクラッチで信頼区間の式を作成するという手間は、それほど無いのではないか。

set.seed(20160322)
my.samples <- 200L  # サンプル数
my.trials <- 100L  # 実験回数
my.conf.p <- 0.95  # 信頼水準
my.mat <- matrix(rnorm(my.samples * my.trials, mean = 0, sd = 1),
                 ncol = my.samples, byrow = TRUE)
my.means <- apply(my.mat, 1, mean)
my.sds <- apply(my.mat, 1, sd)
my.qt <- qt(p = 1 - (1 - my.conf.p) / 2,
            df = my.trials - 1)
my.conf.l <- my.means - my.qt * my.sds / sqrt(my.samples)
my.conf.u <- my.means + my.qt * my.sds / sqrt(my.samples)
plot(x = 1:my.trials, 
     y = my.means,
     xlim = c(0, my.trials + 1),
     ylim = c((min(my.conf.l) - abs(min(my.conf.l)) * 0.05),
              (max(my.conf.u) + abs(min(my.conf.l)) * 0.05)),
     main = paste(sprintf("%d", as.integer(my.conf.p * 100)),
                  "% confidence intervals of ",
                  sprintf("%d", my.trials),
                  " trials from N(0, 1) with ",
                  sprintf("%d", my.samples),
                  " samples", sep = ""),
     xlab = "trials",
     ylab = "means of each trials and its confidence intervals",
     pch = 16,
     cex = 0.5)
segments(x0 = 1:my.trials,
         y0 = my.conf.l,
         y1 = my.conf.u)
abline(h = 0, col = "red")
length(which(my.conf.l < 0 & my.conf.u > 0)) / my.trials
# [1] 0.97

　上記は、エレガントなコードとは言えないし、出力も綺麗ではないなど、コーディング技術上、色々と参考にしたくないと思われてしまうに違いないコードだが、ブログの更新停止状態を打開するためにも、恥を忍んで公表することにした。泥縄式コーディングでも、目的達成の上では、何もないよりも遙かに有効であることを示したかったという意図もある。なお、自分が設定した変数を全部「my」クラス風にしており、見た目に気持ち悪いのは、『サクラエディタ』の「正規表現キーワード機能」でハイライトしているためである。言い訳になるが、「my.」表記は、『R』の多数パッケージにおける変数命名規則が多岐にわたるという事情を斟酌したものである。