頻度論による確率の定義

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


確率分布

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

正規分布 (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のヒストグラム")))

得られたヒストグラムを、上で描いたPDF と比べてみよう。

正規分布以外の確率分布

正規分布に以外の確率分布も同様に扱うことができる。 様々な分布（たとえば、一様分布、ベータ分布、ベルヌーイ分布、二項分布、ポアソン分布など）について、CDF とPDF を図示し、それぞれの分布のパラメタを理解しよう（提出不要の宿題）。


誕生日問題

誕生日問題をR で分析してみよう。Rには、あらかじめ pbirthday() と qbirthday() という関数が用意されているので、これらの関数を使う（詳しい使い方はヘルプを参照）。

まず、集団を構成する人数（横軸）と、その集団内で少なくとも1組が同じ誕生日になる確率（縦軸）の関係を図示してみよう。

## plot the prob of matching birthday vs. n of people
n.vec <- 1:365
probs <- sapply(n.vec, pbirthday)
plot(n.vec, probs, type = "l", col = "tomato", lwd = 2, xlim = c(0, 100),
     xlab = "グループの人数",
     ylab = "確率",
     main = "同じ誕生日のペアが少なくとも1組いる確率：pbirthday()")
abline(h = c(0, 1), col = "gray")

次に、集団内で少なくとも1組が同じ誕生日になる確率をある値（横軸）にするために必要な集団の人数（縦軸）を図示してみよう。

prob.vec <- seq(0, 1, by = 0.01)
n.people <- sapply(prob.vec, qbirthday)
plot(prob.vec, n.people, typ = "l", col = "royalblue", lwd = 2,
     xlab = "少なくとも1組は同じ誕生日になる確率",
     ylab = "グループの人数",
     main = "ある確率を達成するのに必要な人数：qbirthday()")
abline(v = c(0, 1), col = "gray")

このように、直感的な予測よりはだいぶ少ない人数で高い確率が得られることがわかる。



授業の内容に戻る