頻度論による確率の定義

授業用のGitHubリポジトリ を開き、R/frequentisProb.R を実行しなさい。また、自分でシミュレーションを行い、様々な「長期的な相対頻度としての確率」を確認しなさい。


確率分布

Rには主な確率分布を扱うための関数が用意されている。 対象となる確率分布の名前をxxxとすると、 確率密度にはdxxx()(density のd)、累積分布にはpxxx()(probability のp)、累積分布の逆関数には qxxx()(quantile のq)、無作為抽出にはrxxx()(random のr)という関数を使うことができる。 たとえば、正規分布 (normal distribution) を扱うときは、xxxnorm を入れる。

正規分布 (normal distribution) の場合

試しに、正規分布の累積分布関数 (CDF) を図示してみよう。正規分布のCDFは、pnorm() で得られる。 正規分布のパラメタは、平均 \(\mu\) と分散 \(\sigma^2\) であるが、Rでは平均 mean と標準偏差 sd をパラメタとして利用する。デフォルトはそれぞれ0と1(すなわち、標準正規分布)である。また、確率変数 \(X\) の CDFは、\(F_X(x) = \mathrm{Pr}(X \leq x)\) であるので、\(x\) を指定する必要がある。curve() を使って分布曲線を図示すると以下のようになる。

curve(pnorm(x, mean= 50, sd = 10), from = 0, to = 100,
      lwd = 2, col = "royalblue",
      ylab = expression(paste("Pr(",X <= x, ")")),
      main = expression(paste("N(50, ", 10^2, ") の CDF")))
abline(h = c(0, 1), col = "gray")
abline(v = 50, lty = "dashed")
axis(1, 50)

同様に、標準正規分布の確率密度曲線 (PDF) を図示してみよう。標準積分布の PDF は、dnorm() で求める。

curve(dnorm(x, mean = 50, sd = 10), from = 0, to = 100,
      lwd = 2, col = "tomato",
      ylab = "確率密度: p(x)",
      main = expression(paste("N(50, ", 10^2, ") の PDF")))
abline(v = 50, lty = "dashed")
axis(1, 50)

次に、ある \(p\) について、\(F_X(x) = \mathrm{Pr}(X \leq x) = p\) となる \(x\) を与える関数、すなわちCDF の逆関数を使い、\(p=0.025\) となる \(x\)\(p=0.975\) となる\(x\) を同時に求めてみよう。使う関数は qnorm() であり、正規分布のパラメタの他に、確率 \(p\) をベクトルで与える。 求める値は、分布の中心95%をとったときの下限値と上限値である。

qnorm(c(.025, .975), mean = 50, sd = 10)
## [1] 30.40036 69.59964

よって、N(\(50, 10^2\)) の95%は、30.4から69.6 の間にあることがわかる(これは、95%信頼区間ではない。小島『完全独習統計学入門』に登場する言葉を使えば、「95%予言的中区間」である)。

最後に、\(X \sim \mathrm{N}(50, 10^2)\) から1000個の値を無作為に抽出し、それをヒストグラム(縦軸は確率密度)にしてみよう。正規分布の無策抽出に使う関数は、rnorm() である。正規分布のパラメタの他に、抽出する値の個数を指定する。(無策抽出を行う際は、後で結果が再現可能なように、乱数の種 [seed] をset.seed() で設定したほうがよい。Seed にする数はデタラメに 決めればよいが、1つの方法はこの例のように日付を使うことである。)

set.seed(2015-04-15)
x <- rnorm(1000, mean = 50, sd = 10)
hist(x, freq = FALSE, ylab = "確率密度",
     main = expression(paste("N(50, ", 10^2, ") から無作為抽出したxのヒストグラム")))