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))}統計データ解析I
第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
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
1 14 47 151 411 847 1379 1850 1996 1624 1016 469 151 42 2 my_data |> table() |> prop.table() # 確率値に変換my_data
2 3 4 5 6 7 8 9 10 11 12
0.0001 0.0014 0.0047 0.0151 0.0411 0.0847 0.1379 0.1850 0.1996 0.1624 0.1016
13 14 15 16
0.0469 0.0151 0.0042 0.0002 tibble形式のデータフレームを作成するには関数 dplyr::count() を利用して例えば以下のようにすればよい．
my_data |>
as_tibble_col() |> # ベクトルから1列のtibbleを作成．列名の既定値は value
count(value) # value 列の頻度を集計する# A tibble: 15 × 2
value n
<int> <int>
1 2 1
2 3 14
3 4 47
4 5 151
5 6 411
6 7 847
7 8 1379
8 9 1850
9 10 1996
10 11 1624
11 12 1016
12 13 469
13 14 151
14 15 42
15 16 2棒グラフを利用して頻度表を視覚化する．
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: 15 × 3
値 観測値 理論値
<int> <dbl> <dbl>
1 2 0.0001 0.000116
2 3 0.0014 0.000812
3 4 0.0047 0.00396
4 5 0.0151 0.0142
5 6 0.0411 0.0392
6 7 0.0847 0.0840
7 8 0.138 0.142
8 9 0.185 0.189
9 10 0.200 0.198
10 11 0.162 0.162
11 12 0.102 0.101
12 13 0.0469 0.0468
13 14 0.0151 0.0150
14 15 0.0042 0.00301
15 16 0.0002 0.000282package::janitor を利用する場合．
library(janitor)
my_data |> as_tibble_col() |>
tabyl(value) |> # 確率
set_names(c("値","頻度","観測値")) |> # 列名を変更
mutate(理論値 = dbinom(値, size = n, prob = p)) # 理論値を追加 値 頻度 観測値 理論値
2 1 0.0001 0.0001159641
3 14 0.0014 0.0008117488
4 47 0.0047 0.0039572755
5 151 0.0151 0.0142461918
6 411 0.0411 0.0391770274
7 847 0.0847 0.0839507729
8 1379 0.1379 0.1416669293
9 1850 0.1850 0.1888892391
10 1996 0.1996 0.1983337011
11 1624 0.1624 0.1622730282
12 1016 0.1016 0.1014206426
13 469 0.0469 0.0468095274
14 151 0.0151 0.0150459195
15 42 0.0042 0.0030091839
16 2 0.0002 0.0002821110Poisson 分布
問題
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) # 表示領域を指定する(表示しない部分について警告が出る)