library(conflicted) # 関数名の衝突を警告
conflicts_prefer( # 優先的に使う関数を指定
dplyr::filter(),
dplyr::select(),
dplyr::lag(),
)
library(tidyverse)
#' 日本語を用いるので macOS ではフォントの設定を行う
if(Sys.info()["sysname"] == "Darwin") { # macOS か調べて日本語フォントを指定
theme_update(text = element_text(family = "HiraginoSans-W4"))
update_geom_defaults("text", list(family = theme_get()$text$family))
update_geom_defaults("label", list(family = theme_get()$text$family))}第8講 確率分布
いろいろな離散分布と連続分布
準備
以下で利用する共通パッケージを読み込む．
二項分布
問題
二項分布に関して以下を考察せよ．
平均と分散の理論計算を確認しなさい．
X_{1},\dotsc,X_n を成功確率 p の Bernoulli 分布に従う独立同分布な確率変数列とする． このとき \sum_{i=1}^nX_i の分布は，試行回数 n , 成功確率 p の二項分布に従う．これをグラフに描画して確認しなさい．
解答欄
解答例
試行を模擬する関数を定義する．
n <- 16
p <- 0.6
my_random <- function(){ # Bernolli分布をm個生成して合計
sum(rbinom(n, size=1, prob=p))}確率シミュレーションを行う．
mc <- 10000 # 実験回数を指定
my_data <- replicate(mc, my_random())出現値ごとの度数分布表，および正規化した確率表を表示する．
my_data |> table() # base::table による頻度表の作成my_data
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
1 1 10 32 147 406 813 1352 1771 2022 1737 1076 448 152 31 1 my_data |> table() |> prop.table() # 確率値に変換my_data
1 2 3 4 5 6 7 8 9 10 11
0.0001 0.0001 0.0010 0.0032 0.0147 0.0406 0.0813 0.1352 0.1771 0.2022 0.1737
12 13 14 15 16
0.1076 0.0448 0.0152 0.0031 0.0001 tibble形式のデータフレームを作成するには関数 dplyr::count() を利用して例えば以下のようにすればよい．
my_data |>
as_tibble_col() |> # ベクトルから1列のtibbleを作成．列名の既定値は value
count(value) # value 列の頻度を集計する# A tibble: 16 × 2
value n
<int> <int>
1 1 1
2 2 1
3 3 10
4 4 32
5 5 147
6 6 406
7 7 813
8 8 1352
9 9 1771
10 10 2022
11 11 1737
12 12 1076
13 13 448
14 14 152
15 15 31
16 16 1棒グラフを利用して頻度表を視覚化する．
my_data |> as_tibble_col() |>
count(value) |>
ggplot(aes(x = value, y = n)) +
geom_bar(stat = "identity") # 最も単純な表示
理論値(確率)は関数 dbinom() で計算できるので，組み合わせて横並びに表示する．
my_data |> as_tibble_col() |>
count(value) |>
set_names(c("値","観測値")) |> # 列名を変更
mutate(観測値 = 観測値/sum(観測値), # 確率に変換
理論値 = dbinom(値, size = n, prob = p)) |> # 理論値を追加
pivot_longer(!値, values_to = "確率") |>
ggplot(aes(x = 値, y = 確率, fill = name)) +
geom_bar(stat = "identity",
width = 0.8, # 並びがわかりやすいように幅を調整
position = "dodge") +
labs(fill = NULL, # 凡例のfillの名称(name)を消去
title = paste0("二項分布(試行回数", n, ", 成功確率", p, ")"))
Note
上記の例では dplyr::count を利用する場合を示したが， 視覚化のためのデータフレームはいろいろな方法で作成できる．
base::table() および balse::prop.table() を利用する場合．
my_table <- my_data |> table() |> prop.table()
tibble(値 = as.integer(names(my_table)), # 列名は文字列なので整数値に変換(as.numericでも良い)
観測値 = as.vector(my_table), # table型をvector型に変換(as.numericでも良い)
理論値 = dbinom(値, size = n, prob = p))# A tibble: 16 × 3
値 観測値 理論値
<int> <dbl> <dbl>
1 1 0.0001 0.0000103
2 2 0.0001 0.000116
3 3 0.001 0.000812
4 4 0.0032 0.00396
5 5 0.0147 0.0142
6 6 0.0406 0.0392
7 7 0.0813 0.0840
8 8 0.135 0.142
9 9 0.177 0.189
10 10 0.202 0.198
11 11 0.174 0.162
12 12 0.108 0.101
13 13 0.0448 0.0468
14 14 0.0152 0.0150
15 15 0.0031 0.00301
16 16 0.0001 0.000282 package::janitor を利用する場合．
library(janitor)
my_data |> as_tibble_col() |>
tabyl(value) |> # 確率
set_names(c("値","頻度","観測値")) |> # 列名を変更
mutate(理論値 = dbinom(値, size = n, prob = p)) # 理論値を追加 値 頻度 観測値 理論値
1 1 0.0001 1.030792e-05
2 1 0.0001 1.159641e-04
3 10 0.0010 8.117488e-04
4 32 0.0032 3.957275e-03
5 147 0.0147 1.424619e-02
6 406 0.0406 3.917703e-02
7 813 0.0813 8.395077e-02
8 1352 0.1352 1.416669e-01
9 1771 0.1771 1.888892e-01
10 2022 0.2022 1.983337e-01
11 1737 0.1737 1.622730e-01
12 1076 0.1076 1.014206e-01
13 448 0.0448 4.680953e-02
14 152 0.0152 1.504592e-02
15 31 0.0031 3.009184e-03
16 1 0.0001 2.821110e-04Poisson 分布
問題
Poisson 分布に関して以下を考察せよ．
平均と分散の理論計算を確認しなさい．
X,Y を独立な2つの確率変数とし，それぞれ強度 \lambda_{1},\lambda_{2} の Poisson 分布に従うとする． このとき, 和 X+Y の分布は強度 \lambda_{1}+\lambda_{2} の Poisson 分布に従う． これをグラフに描画して確認しなさい．
解答欄
解答例
試行を模擬する関数を定義する．
lambda1 <- 5
lambda2 <- 12
my_random <- function(){ # 2つの Poisson 分布の和
rpois(1, lambda=lambda1)+rpois(1, lambda=lambda2)}確率シミュレーションを行ない，視覚化する．
mc <- 10000
my_data <- replicate(mc, my_random())
my_data |> as_tibble_col() |>
count(value) |> # 頻度を計算
set_names(c("値","観測値")) |> # 列名を変更
mutate(観測値 = 観測値/sum(観測値), # 確率に変換
理論値 = dpois(値, lambda=lambda1+lambda2)) |> # 理論値を追加
pivot_longer(!値, values_to = "確率") |>
ggplot(aes(x = 値, y = 確率, fill = name)) +
geom_bar(stat = "identity",
width = 0.8,
position = "dodge") +
labs(fill = NULL,
title = paste0("Poisson 分布(強度", lambda1+lambda2, ")"))
正規分布
問題
正規分布に関して以下を考察せよ．
U_{1},U_{2} を (0,1) 上の一様分布に従う独立な確率変数とする． このとき \begin{align} X_{1}&=\sqrt{-2\log(U_{1})}\cos(2\pi U_{2}),\\ X_{2}&=\sqrt{-2\log(U_{1})}\sin(2\pi U_{2}) \end{align} とおくと， X_{1},X_{2} は独立かつともに標準正規分布に従う． これをグラフに描画して確認しなさい．
Noteこの変換を Box-Muller 変換と呼ぶ．
解答欄
解答例
試行を模擬する関数を作成する．
my_random <- function(){ # 一方の分布を確認する
u1 <- runif(1)
u2 <- runif(1)
return(sqrt(-2*log(u1))*cos(2*pi*u2))}確率シミュレーションを行い，ヒストグラムを用いて視覚化する．
mc <- 10000 # 実験回数を指定
replicate(mc, my_random()) |> # 確率シミュレーション
as_tibble_col() |> # 既定値では列名は value
ggplot(aes(x = value)) +
geom_histogram(bins = 30, # 既定値では頻度で表示される
fill ="blue",
colour = "white")
理論値は関数 stats::dnorm() で計算できるので， 関数 ggplot2::geom_funciton() を用いて重ね描きする．
replicate(mc, my_random()) |>
as_tibble_col() |>
ggplot(aes(x = value)) +
geom_histogram(aes(y = after_stat(density)), # 確率で表示
bins = 30,
fill ="lightblue",
colour = "white") +
geom_function(fun = dnorm,
colour = "red") +
labs(title = "標準正規分布")
Box-Muller法で作られる2つの確率変数の関係を調べてみる． それぞれを取り出せるように関数を作成する．
boxmuller <- function(){
u1 <- runif(1)
u2 <- runif(1)
x1 <- sqrt(-2*log(u1))*cos(2*pi*u2)
x2 <- sqrt(-2*log(u1))*sin(2*pi*u2)
return(c(x1,x2))
}2つの変数の散布図を描く．
x <- replicate(mc, boxmuller()) # 2行xmc列の行列が得られる
tibble(x1 = x[1,], x2 = x[2,]) |>
ggplot(aes(x = x1, y = x2)) +
geom_point(colour = "royalblue")
x1,x2 は同じ分布に従う独立な変数なので以下ではまとめて扱う． 個別に扱う場合は x を x[1,] などとすれば良い．
mu <- mean(x)
sigma <- sd(x)
tibble(x = as.vector(x)) |>
ggplot(aes(x = x)) +
geom_histogram(aes(y = after_stat(density)),
bins = 30,
fill ="lightblue",
colour = "white") +
geom_function(fun = dnorm,
args = list(mean = mu, sd = sigma),
colour = "red") +
labs(title = paste0("正規分布(平均", round(mu,2), # 2桁で四捨五入
", 分散", round(sigma^2,2), ")"))
Note
散布図を描くためのデータフレームは以下のようにして作ることもできる． 高次元の場合には関数を利用した方が簡単に記述できる．
x |> t() |>
as_tibble(.name_repair = "minimal") |>
set_names(paste0("x", 1:2)) |>
ggplot(aes(x = x1, y = x2)) +
geom_point(colour = "royalblue")
\chi^2 分布
問題
\chi^2 分布に関して以下を考察せよ．
- 標準正規分布に従う k 個の独立な確率変数の二乗和は自由度 k の \chi^2 分布に従うことを確認しなさい．
解答欄
解答例
試行を模擬する関数を定義する．
k <- 8 # 自由度を設定
my_random <- function(){
sum(rnorm(k)^2)} # k個の標準正規乱数の二乗和確率シミュレーションを行い，視覚化する．
mc <- 10000 # 実験回数を指定
my_data <- replicate(mc, my_random()) # 注で使うので保存しておく
my_data |>
as_tibble_col() |>
ggplot(aes(x = value)) +
geom_histogram(aes(y = after_stat(density)),
bins = 30,
fill ="lightblue",
colour = "white") +
geom_function(fun = dchisq,
args = list(df = k),
colour = "red") +
labs(title = bquote(paste(chi^2,"分布 (自由度",.(k),")")))
Note
関数 bquote() は関数 expression() と同様な働きをするが，自由度の k を .(k) で評価して取り込むことができる． 以下と表示を比較してみよう．
my_data |>
as_tibble_col() |>
ggplot(aes(x = value)) +
geom_histogram(aes(y = after_stat(density)),
bins = 30,
fill ="lightblue",
colour = "white") +
geom_function(fun = dchisq,
args = list(df = k),
colour = "red") +
labs(title = expression(paste(chi^2,"分布 (自由度",k,")"))) # この行だけ違う
t 分布
問題
t 分布に関して以下を考察せよ．
- Z を標準正規分布に従う確率変数，Y を自由度 k の \chi^2 分布に従う確率変数とし，Z,Y は独立であるとする． このとき確率変数 \begin{equation} \frac{Z}{\sqrt{Y/k}} \end{equation} は自由度 k の t 分布に従うことを確認しなさい．
解答欄
解答例
試行を模擬する関数を作成する．
k <- 8 # 自由度を設定
my_random <- function(){
sum(rnorm(k)^2)} # k個の標準正規乱数の二乗和確率シミュレーションを行い，視覚化する．
mc <- 10000 # 実験回数を指定
my_data <- replicate(mc, my_random()) # 注のところでもう一度描くので保存しておく
my_data |>
as_tibble_col() |>
ggplot(aes(x = value)) +
geom_histogram(aes(y = after_stat(density)),
bins = 30,
fill ="lightblue",
colour = "white") +
geom_function(fun = dchisq,
args = list(df = k),
colour = "red") +
labs(title = bquote(paste(chi^2,"分布 (自由度",.(k),")")))
Note
関数 bquote() は関数 expression() と同様な働きをするが， 自由度の k を .(k) で評価して取り込むことができる． 以下と比較してみよう．
my_data |>
as_tibble_col() |>
ggplot(aes(x = value)) +
geom_histogram(aes(y = after_stat(density)),
bins = 30,
fill ="lightblue",
colour = "white") +
geom_function(fun = dchisq,
args = list(df = k),
colour = "red") +
labs(title = expression(paste(chi^2,"分布 (自由度",k,")")))
F 分布
問題
F 分布に関して以下を考察せよ．
- Y_{1},Y_{2} をそれぞれ自由度 k_{1},k_{2} の \chi^2 分布に従う独立な確率変数とする． このとき, 確率変数 \begin{equation} \frac{Y_{1}/k_{1}}{Y_{2}/k_{2}} \end{equation} は自由度 k_{1},k_{2} の F 分布に従うことを確認しなさい．
解答欄
解答例
試行を模擬する関数を作成する．
k1 <- 20
k2 <- 10
my_random <- function(){
y1 <- rchisq(1, df=k1) # 自由度20のカイ二乗分布
y2 <- rchisq(1, df=k2) # 自由度10のカイ二乗分布
## y1 <- sum(rnorm(k1)^2) # 正規乱数を用いてもよい
## y2 <- sum(rnorm(k2)^2)
return((y1/k1)/(y2/k2))}確率シミューレションを行い，視覚化する．
mc <- 10000 # 実験回数を指定
replicate(mc, my_random()) |>
as_tibble_col() |>
ggplot(aes(x = value)) +
geom_histogram(aes(y = after_stat(density)),
bins = 30,
fill ="lightblue",
colour = "white") +
geom_function(fun = df,
args = list(df1 = k1, df2 = k2),
colour = "red") +
labs(title = bquote(paste(frac(Y[1]/k[1],Y[2]/k[2]),
" (",k[1]==.(k1),
", ",k[2]==.(k2),")")))
Note
グラフの一部に着目したい場合は xlim/ylim で表示を調整することができる．
last_plot() + # 直前のプロット
xlim(0,6) + ylim(0,0.9) # 表示領域を指定する(表示しない部分について警告が出る)