====== 【社会情報処理】第3回 様々な確率分布 ======
統計学では母集団から標本を抽出し、その標本から母集団についての推定や検定を行う。
そのためには確率分布が必要となる。

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

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

{{ socialip_2021:distribution.png?nolink |正規分布の例}}

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

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

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

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

^  関数                                              ^  機能                                                                                         ^
| ''dunif(x, min = 0, max = 1)''                     | 確率密度関数 $f(x)$ である。                                                                  |
| ''punif(q, min = 0, max = 1, lower.tail = TRUE)''  | 下側確率 $P(X\le q)$ である。\\ 引数 ''lower.tail = FALSE'' なら上側確率になる。              |
| ''qunif(p, min = 0, max = 1, lower.tail = TRUE)''  | 下側確率が ''p'' となるパーセント点である。\\ 引数 ''lower.tail=FALSE'' なら上側確率になる。  |
| ''runif(n, min = 0, max = 1)''                     | 一様分布の乱数 ''n'' 個を生成する。\\ 乱数の範囲は $\mathrm{min}<x<\mathrm{max}$ である。     |
引数の中に「''=''」とあるのは、その引数を省略した場合のデフォルト値を表す。


例えば $1<x<7$ の一様分布の乱数（__**一様乱数**__）は以下で発生させる。
<code rsplus>
> 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
</code>
整数の一様乱数が欲しければ、''as.integer()'' と組み合わせて小数点以下を切り捨てる。

1 から 6 までの一様乱数（サイコロの目）が欲しければ以下のようにする。
<code rsplus>
> as.integer(runif(10, 1, 7))
 [1] 4 1 5 3 6 6 5 5 6 4
</code>

^  一様乱数 $x$  ^  ''as.integer(x)''  ^
|  $1<x<2$       |  $1$                |
|  $2\le x<3$    |  $2$                |
|  $3\le x<4$    |  $3$                |
|  $4\le x<5$    |  $4$                |
|  $5\le x<6$    |  $5$                |
|  $6\le x<7$    |  $6$                |

''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'' で指定する。

^  関数                                    ^  機能                                                                                         ^
| ''dbinom(x, size, prob)''                  | 確率密度関数 $f(x)$ である。                                                                  |
| ''pbinom(q, size, prob, lower.tail = TRUE)''  | 下側確率 $P(X\le q)$ である。\\ 引数 ''lower.tail = FALSE'' なら上側確率になる。                |
| ''qbinom(p, size, prob, lower.tail = TRUE)''  | 下側確率が ''p'' となるパーセント点である。\\ 引数 ''lower.tail = FALSE'' なら上側確率になる。  |
| ''rbinom(n, size, prob)''                  | 二項分布の乱数 ''n'' 個を生成する。                                               |


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

<file rsplus dist_binomial.R[enable_line_numbers="true"]>
# 二項分布

# パラメータ
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)
</file>

''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'' で指定する。

^  関数                                     ^  機能                                                                                         ^
| ''dnorm(x, mean = 0, sd = 1)''                  | 確率密度関数 $f(x)$ である。                                                                  |
| ''pnorm(q, mean = 0, sd = 1, lower.tail = TRUE)''  | 下側確率 $P(X\le q)$ である。\\ 引数 ''lower.tail = FALSE'' なら上側確率になる。                |
| ''qnorm(p, mean = 0, sd = 1,lower.tail = TRUE)''  | 下側確率が ''p'' となるパーセント点である。\\ 引数 ''lower.tail = FALSE'' なら上側確率になる。  |
| ''rnorm(n, mean = 0, sd = 1)''                  | 正規分布の乱数 ''n'' 個を生成する。                                               |

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

<file rsplus dist_normal.R[enable_line_numbers="true"]>
# 正規分布

# パラメータ
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)
</file>

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

===== 小テスト =====

[[KMS>|Moodle Server（非公式）]]で第3回の小テストを受験しなさい。
===== スチューデントのt分布 =====
自由度 $\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'' で指定する。

^  関数                             ^  機能                                                                                           ^
| ''dt(x, df)''                     | 確率密度関数 $f(x)$ である。                                                                    |
| ''pt(q, df, lower.tail = TRUE)''  | 下側確率 $P(X\le q)$ である。\\ 引数 ''lower.tail = FALSE'' なら上側確率になる。                |
| ''qt(p, df, lower.tail = TRUE)''  | 下側確率が ''p'' となるパーセント点である。\\ 引数 ''lower.tail = FALSE'' なら上側確率になる。  |
| ''rt(n, df)''                     | スチューデントの t 分布の乱数 ''n'' 個を生成する。                                              |


次のプログラム dist_t.R は $t(\nu)$ のグラフを描く。
<file rsplus dist_t.R[enable_line_numbers="true"]>
# スチューデントの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)
</file>

===== カイ二乗分布 =====
自由度 $\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'' で指定する。

^  関数                             ^  機能                                                                                         ^
| ''dchisq(x, df)''                  | 確率密度関数 $f(x)$ である。                                                                  |
| ''pchisq(q, df, lower.tail = TRUE)''  | 下側確率 $P(X\le q)$ である。\\ 引数 ''lower.tail = FALSE'' なら上側確率になる。                |
| ''qchisq(p, df, lower.tail = TRUE)''  | 下側確率が ''p'' となるパーセント点である。\\ 引数 ''lower.tail = FALSE'' なら上側確率になる。  |
| ''rchisq(n, df)''                  | カイ二乗分布の乱数 ''n'' 個を生成する。                                           |

次のプログラム dist_chisq.R は $\chi^{2}(\nu)$ のグラフを描く。
<file rsplus dist_chisq.R[enable_line_numbers="true"]>
# カイ二乗分布

# パラメータ
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)
</file>

===== F分布 =====
自由度 $(\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'' でそれぞれ指定する。

^  関数                             ^  機能                                                                                         ^
| ''df(x, df1, df2)''                  | 確率密度関数 $f(x)$ である。                                                                  |
| ''pf(q, df1, df2, lower.tail = TRUE)''  | 下側確率 $P(X\le q)$ である。\\ 引数 ''lower.tail = FALSE'' なら上側確率になる。                |
| ''qf(p, df1, df2, lower.tail = TRUE)''  | 下側確率が ''p'' となるパーセント点である。\\ 引数 ''lower.tail = FALSE'' なら上側確率になる。  |
| ''rf(n, df1, df2)''                  | F分布の乱数 ''n'' 個を生成する。                                           |

次のプログラム dist_f.R は $F(\nu_{1},\nu_{2})$ のグラフを描く。
<file rsplus dist_f.R[enable_line_numbers="true"]>
# 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)
</file>
===== 標本抽出（サンプリング） =====
母集団から無作為の__**標本抽出（サンプリング）**__には一様乱数を使う。
一様乱数ならば ''runif()'' があるが、標本抽出を行うための関数 ''sample()'' があるのでこれを利用すればよい。

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

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

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

<file rsplus sampling.R[enable_line_numbers="true"]>
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)
</file>



===== 課題  =====
[[KMS>|Moodle Server（非公式）]] で第3回の課題を行いなさい。

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

----
==== 解答例 ====