統計学では母集団から標本を抽出し、その標本から母集団についての推定や検定を行う。 そのためには確率分布が必要となる。

R では様々な確率分布を扱うための関数が用意されている。 関数は頭文字によって機能が分けられ、頭文字が d の関数は確率密度関数、p の関数は累積確率、q の関数はパーセント点、r の関数はその確率分布の乱数を返す。

例えば、正規分布の場合の関数との対応は以下の図のようになる。

ここでは代表的な確率分布をいくつか挙げる。

一様分布は $x$ の区間 $[a,b]$ で確率が一定の確率分布である。

確率密度関数 \[ f(x)=\frac{1}{b-a} \]

一様分布に関する R の関数は以下の通り。 $a$ は引数 min、$b$ は引数 max で指定する。

引数の中に「=」とあるのは、その引数を省略した場合のデフォルト値を表す。

例えば $1<x<7$ の一様分布の乱数（一様乱数）は以下で発生させる。

> runif(10, 1, 7)
 [1] 4.580631 1.360485 5.536313 3.697583 6.004030 6.948436 5.668975 5.357032 6.656362
[10] 4.427038

整数の一様乱数が欲しければ、as.integer() と組み合わせて小数点以下を切り捨てる。

1 から 6 までの一様乱数（サイコロの目）が欲しければ以下のようにする。

> as.integer(runif(10, 1, 7))
 [1] 4 1 5 3 6 6 5 5 6 4

runif() で最小値の値は乱数として発生しないので $x=1$ を含まないが、ほとんど違いはないので気にしなくていい。

2種類の可能な結果（「成功」と「失敗」とする）を生じる実験、あるいは観測があり、これを同じ条件で繰り返すことをベルヌーイ試行という。 ベルヌーイ試行を $1$ 回行って「成功」が出る確率を $p$ としたとき、そのベルヌーイ試行を $n$ 回くり返して「成功」が $x$ 回起こる確率の確率分布が二項分布である。 二項分布は $Bi(n,p)$ で表す。

確率密度関数 \[ f(x)=_{n}C_{x}p^{x}(1-p)^{n-x} \] ここで $_{n}C_{x}$ は二項係数である。 \[ _{n}C_{x}=\frac{n!}{x!(n-x)!} \] $Bi(n,p)$ は離散型の確率分布であり、$x$ は $0,1,2,\cdots, n$ の値しかとらない。

$Bi(n,p)$ に関する R の関数は以下の通り。 $n$ は引数 size、$p$ は引数 prob で指定する。

次のプログラム dist_binomial.R は $Bi(n,p)$ のグラフを描き、確率を出力する。

dist_binomial.R

# 二項分布
 
# パラメータ
n <- 10
p <- 0.5
 
# 確率密度関数
x <- 0 : n
f <- dbinom(x, n, p)
 
# グラフ
title <- sprintf("Bi(%d,%g)", n, p)
barplot(f, names.arg = x, xlab = "x", ylab = "f", space = 0, main = title)
 
# 値の出力
print(sprintf("f(%d)=%f", x, f), quote = FALSE)

barplot() は棒グラフを描く関数である。

平均 $\mu$、分散 $\sigma^{2}$ の正規分布を $N(\mu,\sigma^{2})$ で表す。 特に $\mu=0$、$\sigma^{2}=1$ の正規分布 $N(0,1)$ を標準正規分布 という。

確率密度関数 \[ f(x)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left[-\frac{(x-\mu)^{2}}{2\sigma^{2}}\right] \]

$N(\mu,\sigma^{2})$ に関する R の関数は以下の通り。 平均 $\mu$ は引数 mean、標準偏差 $\sigma$ は引数 sd で指定する。

次のプログラム dist_normal.R は $N(\mu,\sigma^{2})$ のグラフを描く。

dist_normal.R

# 正規分布
 
# パラメータ
mu <- 5
sigma <- 10
 
# 確率密度関数のグラフ
xmin <- mu - sigma * 5
xmax <- mu + sigma * 5
title <- sprintf("N(%g,%g)", mu, sigma ^ 2)
curve(dnorm(x, mu, sigma), from = xmin, to = xmax, xlab = "x", ylab = "f",
      main = title)

curve() は曲線を描く関数である。

Moodle Server（非公式）で第3回の小テストを受験しなさい。

自由度 $\nu$ のスチューデントのt分布を $t(\nu)$ で表す。

確率密度関数 \[ f(t)=\frac{\Gamma\left(\frac{\nu+1}{2}\right)}{\Gamma\left(\frac{\nu}{2}\right)\sqrt{\nu\pi}}\left(1+\frac{t^{2}}{\nu}\right)^{-\frac{\nu+1}{2}} \] ここで $\Gamma(x)$ はガンマ関数である。

$t(\nu)$ に関する R の関数は以下の通り。 自由度 $\nu$ は引数 df で指定する。

次のプログラム dist_t.R は $t(\nu)$ のグラフを描く。

dist_t.R

# スチューデントのt分布
 
# パラメータ
nu <- 3
 
# 確率密度関数のグラフ
tmin <- -5
tmax <- 5
title <- sprintf("t(%d)", nu)
curve(dt(x, nu), from = tmin, to = tmax, xlab = "t", ylab = "f", main = title)

自由度 $\nu$ のカイ二乗分布を $\chi^{2}(\nu)$ で表す。

確率密度関数 \[ f(\chi^{2})=\frac{1}{2\Gamma\left(\frac{\nu}{2}\right)}\left(\frac{\chi^{2}}{2}\right)^{\frac{\nu}{2}-1}e^{-\frac{\chi^{2}}{2}} \] ここで $\Gamma(x)$ はガンマ関数である。

$\chi^{2}(\nu)$ に関する関数は以下の通り。 自由度 $\nu$ は引数 df で指定する。

次のプログラム dist_chisq.R は $\chi^{2}(\nu)$ のグラフを描く。

dist_chisq.R

# カイ二乗分布
 
# パラメータ
nu <- 5
 
# 確率密度関数のグラフ
cmin <- 0
cmax <- nu * 3
title <- parse(text = sprintf("chi^2*(%d)", nu))
curve(dchisq(x, nu), from = cmin, to = cmax, main = title)

自由度 $(\nu_{1},\nu_{2})$ のF分布を $F(\nu_{1},\nu_{2})$ で表す。

確率密度関数 \[ f(F)=\frac{\nu_{1}^{\frac{\nu_{1}}{2}}\nu_{2}^{\frac{\nu_{2}}{2}}}{B(\frac{\nu_{1}}{2},\frac{\nu_{2}}{2})}F^{\frac{\nu_{1}}{2}-1}(\nu_{2}+\nu_{1}F)^{-\frac{\nu_{1}+\nu_{2}}{2}} \]

$F(\nu_{1},\nu_{2})$ に関する R の関数は以下の通り。 自由度 $(\nu_{1},\nu_{2})$ は引数 df1,df2 でそれぞれ指定する。

次のプログラム dist_f.R は $F(\nu_{1},\nu_{2})$ のグラフを描く。

dist_f.R

# F分布
 
# パラメータ
nu1 <- 4
nu2 <- 5
 
# 確率密度関数のグラフ
fmin <- 0
fmax <- 5
title <- sprintf("F(%d,%d)", nu1, nu2)
curve(df(x, nu1, nu2), from = fmin, to = fmax, main = title)

母集団から無作為の標本抽出（サンプリング）には一様乱数を使う。 一様乱数ならば runif() があるが、標本抽出を行うための関数 sample() があるのでこれを利用すればよい。

> sample(母集団のデータ, n)

母集団のデータはベクトルの必要がある。 n は標本のサイズである。 n を省略すると母集団のすべての要素をランダムに並べて返す。

デフォルトでは非復元抽出の標本抽出になるが、 引数に replace = TRUE を入れると復元抽出になる。

次のプログラム sampling.R は第2回に e-Stat からダウンロードした山形県の市町村データ（ファイル名は yamagata.csv にしている）において、2015年度の人口増減率を母集団として標本抽出を行うプログラムである。 分かりやすいようにデータはソートして出力している。

sampling.R

library(readr)
 
# データの読み込み
df <- read_csv("yamagata.csv", locale = locale(encoding = "SHIFT-JIS"))
df$`/項目` <- NULL
df$調査年 <- as.numeric(sapply(df$調査年, gsub, pattern="年度", replacement = ""))
df <- df[order(df$調査年),]
 
# 2015年度の人口増減率を母集団とする
data <- df[df$調査年 == 2015,][["#A05101_人口増減率【％】"]]
 
print("母集団", quote = FALSE)
print(sort(data))
 
# 標本サイズ
n <- 20
 
print("非復元抽出", quote = FALSE)
sdata1 <- sample(data, n)
print(sort(sdata1), quote = FALSE)
 
print("復元抽出", quote = FALSE)
sdata2 <- sample(data, n, replace = TRUE)
print(sort(sdata2), quote = FALSE)

Moodle Server（非公式） で第3回の課題を行いなさい。

締め切り：2021年7月8日（木）20時