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
第10講 練習問題 解答例
準備
以下で利用する共通パッケージを読み込む．
一様分布の平均の推定
問題
X を一様乱数に従う確率変数とし，平均値の推定量として以下を考える． それぞれの推定量の分散を比較しなさい．
- 標本平均 (mean)
- 中央値 (median)
- 最大値と最小値の平均 ((max+min)/2)
ヒント
以下のような関数を作り，確率シミュレーションを行えばよい．
estimate_means <- function(n, min, max){ # 観測データ数
x <- runif(n, min = min, max = max) # 一様乱数を生成，範囲は引数から
return(c(xbar = mean(x), med = median(x), mid = (max(x)+min(x))/2))
} # 3つまとめて計算する関数解答例
平均値の推定を行う関数(標本平均，中央値，最大最小の平均)を定義する．
estimate_means <- function(n, min, max){
x <- runif(n, min = min, max = max)
return(c(xbar = mean(x), # 標本平均
med = median(x), # 中央値
mid = (max(x)+min(x))/2)) # 最大最小の平均
}実験の設定を宣言する．
n <- 10 # 観測データのサンプル数
mc <- 10000 # 実験回数
a <- 0; b <- 50 # 一様乱数の区間確率シミュレーションを行う．
mc_data <-
replicate(mc, estimate_means(n, min = a, max = b)) |> # 3 x mc 型行列
t() |> # 転置して mc x 3 型行列に変換
as_tibble() # データフレーム化
mc_data # 実験結果の一部を確認# A tibble: 10,000 × 3
xbar med mid
<dbl> <dbl> <dbl>
1 31.7 31.3 31.6
2 29.2 29.4 27.3
3 26.5 22.8 23.8
4 28.3 29.2 29.2
5 22.4 24.7 24.2
6 20.5 18.0 25.3
7 20.3 17.9 20.1
8 31.5 29.1 28.7
9 25.9 28.9 25.0
10 27.3 26.5 24.5
# ℹ 9,990 more rowsそれぞれの推定量の平均と分散を計算し，結果を集計する．
mc_data |>
summarize(across(everything(),
list(mean = mean, var = var)))# A tibble: 1 × 6
xbar_mean xbar_var med_mean med_var mid_mean mid_var
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 25.0 21.0 25.0 47.3 25.0 9.58別々のヒストグラムにして比較してみる．
mc_data |>
pivot_longer(everything()) |>
ggplot(aes(x = value)) +
geom_histogram(aes(fill = name),
show.legend = FALSE) + # 凡例を表示しない
facet_wrap(~name) # 横に並べる 
一様分布においては，最大最小の平均が良い推定量であることがわかる．
Note
分位点を使ってもう少し詳しくみてみる．
mc_data |> summary() # 四分位点を表示 (簡便な方法) xbar med mid
Min. :10.22 Min. : 3.413 Min. :11.25
1st Qu.:21.83 1st Qu.:20.068 1st Qu.:23.27
Median :24.90 Median :24.945 Median :25.01
Mean :24.99 Mean :24.971 Mean :25.00
3rd Qu.:28.15 3rd Qu.:29.834 3rd Qu.:26.73
Max. :43.16 Max. :47.231 Max. :39.16 mc_data |> # データフレームとして取得する方法 (以下は一例)
pivot_longer(everything()) |> # long format に変更
group_by(name) |> # 推定量ごとに集計
summarize(as_tibble_row(quantile(value))) # 計算結果を行に並べる# A tibble: 3 × 6
name `0%` `25%` `50%` `75%` `100%`
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 med 3.41 20.1 24.9 29.8 47.2
2 mid 11.3 23.3 25.0 26.7 39.2
3 xbar 10.2 21.8 24.9 28.2 43.2mc_data |> # 別形式のデータフレームの作成例
pivot_longer(everything()) |> # long format に変更
group_by(name) |> # 推定量ごとに集計
reframe(qs=quantile(value), probs=seq(0,1,.25)) # A tibble: 15 × 3
name qs probs
<chr> <dbl> <dbl>
1 med 3.41 0
2 med 20.1 0.25
3 med 24.9 0.5
4 med 29.8 0.75
5 med 47.2 1
6 mid 11.3 0
7 mid 23.3 0.25
8 mid 25.0 0.5
9 mid 26.7 0.75
10 mid 39.2 1
11 xbar 10.2 0
12 xbar 21.8 0.25
13 xbar 24.9 0.5
14 xbar 28.2 0.75
15 xbar 43.2 1 それぞれのヒストグラムを描くこともできる．
for(name in c("xbar","med","mid")) {
name <- sym(name)
gg <-
mc_data |>
ggplot(aes(x = !!name)) +
geom_histogram(bins = 40,
fill = "pink",
colour = "brown") +
lims(x = c(a,b), y = c(0,2500)) + # 同じ大きさの図にする
labs(x = "estimate",
title = paste(name, "の分布"))
print(gg)
}
Note
違う分布で試してみる．
estimate_means2 <- function(n){
x <- rt(n, df = 2) # 自由度2のt分布 (裾が重く平均が推定しにくい分布)
return(c(xbar = mean(x), med = median(x), mid = (max(x)+min(x))/2))
}
mc_data2 <-
replicate(mc, estimate_means2(n)) |> t() |> as_tibble()
#' それぞれの分布を書いてみる
mc_data2 |>
pivot_longer(everything()) |>
ggplot(aes(x = value)) +
geom_histogram(aes(fill = name),
show.legend = FALSE) +
xlim(-8,8) + # x軸を調整する
facet_wrap(~name) 
この分布では中央値が良い推定量となることがわかる．
Tip
Rで用いることのできる色の名前は関数 colors()/colours() で確認することができる. 例えば以下のようにすると視覚的に確認することができる.
cols <- colors()[grep("(grey|[0-1,3-9])", # greyおよび2以外の数を含む色を排除
colors(), invert = TRUE)]
ncols <- 6; nrows <- ceiling(length(cols)/ncols) # 必要なタイルの枚数の計算
tibble(x = rep(1:ncols, nrows),
y = rep(1:nrows, each = ncols),
z = factor(1:(nrows*ncols))) |> # タイルの配置のためのデータフレーム
slice_head(n = length(cols)) |> # 必要な行のみ残す
ggplot(aes(x = x, y = y, fill = z)) +
scale_fill_manual(values = cols) + # fillの色を指定
geom_tile(show.legend = FALSE) + # 色のタイルを配置
geom_text(aes(label = cols), size = 2) # 色名を記入
ガンマ分布による風速データのモデル化
問題
東京都の気候データ (tokyo_weather.csv) の風速 (wind) の項目について以下の問に答えよ．
- 全データを用いてヒストグラム(密度)を作成しなさい．
- ガンマ分布でモデル化して最尤推定を行いなさい．
- 推定した結果をヒストグラムに描き加えて比較しなさい．
Note
ガンマ分布の最尤推定量は以下のようにして作成できる
library("stats4") # 関数mleを利用するため
#' 数値最適化のためには尤度関数を最初に評価する初期値が必要
mle.gamma <- function(x, # 観測データ
nu0 = 1, alpha0 = 1){ # nu, alphaの初期値
## 負の対数尤度関数を定義 (最小化を考えるため)
ll <- function(nu, alpha) # nuとalphaの関数として定義
suppressWarnings(-sum(dgamma(x, nu, alpha, log = TRUE)))
## suppressWarnings は定義域外で評価された際の警告を表示させない
## 最尤推定(負の尤度の最小化)
est <- mle(minuslogl = ll, # 負の対数尤度関数
start = list(nu = nu0, alpha = alpha0), # 初期値
method = "BFGS", # 最適化方法 (選択可能)
nobs = length(x)) # 観測データ数
return(coef(est)) # 推定値のみ返す
}
Note
余力のある者は，自身で収集したデータを用いてモデル化と最尤推定を試みよ．
解答例
データを読み込む．
tw_data <- read_csv("data/tokyo_weather.csv")風速データのヒストグラムを描く．
tw_data |>
ggplot(aes(x = wind)) +
geom_histogram(aes(y = after_stat(density)), # 密度表示
bins = 30,
fill = "skyblue",
colour = "slateblue") +
labs(x = "風速 [m/s]", y = "密度", title = "風速のヒストグラム")
ガンマ分布の最尤推定量を求めるための関数を作成する．
library("stats4") # 関数mleを利用
mle.gamma <- function(x, # 観測データ
nu0 = 1, # nu (shape) の初期値
alpha0 = 1, # alpha (rate) の初期値
verbose = FALSE){ # debug用に追加
## 負の対数尤度関数を定義 (最小化を考えるため)
ll <- function(nu, alpha) # nuとalphaの関数として定義
suppressWarnings(-sum(dgamma(x,
shape = nu,
rate = alpha,
log = TRUE))) # 対数密度
## suppressWarnings は定義域外で評価された際の警告を表示させない
## 最尤推定(負の尤度の最小化)
est <- mle(minuslogl = ll, # 負の対数尤度関数
start = list(nu = nu0, alpha = alpha0), # 初期値
method = "BFGS", # 最適化方法 (選択可能)
nobs = length(x)) # 観測データ数
if(verbose) {
return(est) # verbose=TRUEならmleの結果を全て返す
} else {
return(coef(est)) # 推定値のみ返す
}
}最尤推定を行う．
(theta <- with(tw_data, mle.gamma(wind))) nu alpha
10.919235 4.012087 得られた結果を前出のヒストグラム上に重ね描きする．
last_plot() +
geom_function(fun = dgamma,
args = list(shape = theta[1],
rate = theta[2]),
colour = "orange",
linewidth = 1.2)
Note
シミュレーションにより推定量の良さの検証(一致性や不偏性の確認)を行うと以下のようになる．
set.seed(5678) # 乱数のシード値の指定 (変更して再現性を確認してみよ)
nu <- 5; alpha <- 2 # 真のパラメタ
mc <- 1000 # 実験回数 (計算が重いので少なめにしている)
for(n in c(10, 50, 100)){ # データ数を変えて実験
## Monte-Carlo実験
mc_data <- # 推定値のデータフレーム
replicate(mc, mle.gamma(rgamma(n, nu, alpha))) |>
t() |> as_tibble()
## 結果を密度推定で表示
gg <-
mc_data |>
ggplot(aes(x = nu)) + # nu の推定値の分布
geom_density(fill = "skyblue1",
colour = "skyblue4") +
geom_vline(xintercept = nu, # nu の真値
colour = "tomato",
linewidth = 1.2) +
xlim(0,20) +
labs(x = expression(nu), title = paste("n =",n))
print(gg)
gg <-
mc_data |>
ggplot(aes(x = alpha)) + # alpha の推定値の分布
geom_density(fill = "seagreen1",
colour = "seagreen4") +
geom_vline(xintercept = alpha, # alpha の真値
colour = "tomato",
linewidth = 1.2) +
xlim(0,10) +
labs(x = expression(alpha), title = paste("n =",n))
print(gg)
}
日射量データの区間推定
問題
東京都の気候データ (tokyo_weather.csv) の日射量 (solar) の項目について以下の問に答えよ．
- 全データによる平均値を計算しなさい．
- ランダムに抽出した50点を用いて，平均値の0.9(90%)信頼区間を求めなさい．
- 上記の推定を100回繰り返した際，真の値(全データによる平均値)が信頼区間に何回含まれるか確認しなさい．
Note
余力のある者は，自身で収集したデータで区間推定を試みよ．
解答例
全データによる平均値の計算を行う．
(mu <- tw_data |> pull(solar) |> mean()) # 真値[1] 13.88866これを真値として信頼区間の精度の検証を行う．
ランダムに選んだデータ50点を用いて90%信頼区間を構成する．
set.seed(1357) # 乱数のシード値の指定 (変更して再現性を確認してみよ)
n <- 50
est <- # 平均と標準偏差の推定を行う
tw_data |>
slice_sample(n = n) |> # ランダムにn個取り出す
summarize(across(solar, list(mean = mean, sd = sd)))
xbar <- est[["solar_mean"]] # 標本平均
sigma <- est[["solar_sd"]] # 標本標準偏差
z95 <- qnorm(0.95) # 標準正規分布の0.95分位点
(ci <- c(L = xbar-z95*sigma/sqrt(n),
U = xbar+z95*sigma/sqrt(n))) # 信頼区間 L U
12.66197 15.93843 これを複数回行い，信頼区間の正答率を評価する．
mc <- 100
mc_trial <- function(n){ # nを変えて実験できるように
est <- tw_data |> slice_sample(n = n) |>
summarize(across(solar, list(mean = mean, sd = sd)))
xbar <- est[["solar_mean"]] # 標本平均
sigma <- est[["solar_sd"]] # 標本標準偏差
return(c(L = xbar-z95*sigma/sqrt(n),
U = xbar+z95*sigma/sqrt(n))) # 信頼区間
}
mc_data <- # 信頼区間のMonte-Carlo実験
replicate(mc, mc_trial(n)) |> t() |> as_tibble() |>
mutate(answer = L < mu & mu < U) # 真値が信頼区間に含まれるか
mc_data |> count(answer) # 頻度を見る# A tibble: 2 × 2
answer n
<lgl> <int>
1 FALSE 4
2 TRUE 96
Note
信頼区間について更に多数で評価してみる．
mc <- 2000
mc_data <- # 信頼区間のMonte-Carlo実験
replicate(mc, mc_trial(n)) |> t() |> as_tibble() |>
mutate(answer = L < mu & mu < U)
mc_data |> count(answer) |> mutate(prob=n/sum(n)) # 確率を見る# A tibble: 2 × 3
answer n prob
<lgl> <int> <dbl>
1 FALSE 194 0.097
2 TRUE 1806 0.903
Note
日射量の平均の推定における観測データを全て使った推定値(真値に相当)と信頼区間の関係をグラフを描いて確認してみる．
k <- 20
mc_data |>
slice_sample(n = k) |> # k個ランダムに選んで描く
rowid_to_column(var = "index") |>
ggplot(aes(x = index)) +
geom_errorbar(aes(ymin = L, ymax = U),
colour = "blue",
width = 0.2) +
geom_hline(yintercept = mu, # 観測データを全て使った推定値
colour = "orange",
linewidth = 1.1)