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))}第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() # データフレーム化
print(mc_data) # 実験結果の一部を確認# A tibble: 10,000 × 3
xbar med mid
<dbl> <dbl> <dbl>
1 18.8 16.2 24.2
2 17.7 11.5 24.7
3 23.1 23.5 25.5
4 24.3 23.2 26.0
5 19.1 14.8 25.0
6 24.9 29.8 21.8
7 27.1 28.8 26.7
8 22.8 27.7 23.0
9 24.2 23.6 25.7
10 33.9 36.1 31.8
# ℹ 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.6 25.0 9.36別々のヒストグラムにして比較してみる．
mc_data |>
pivot_longer(everything()) |>
ggplot(aes(x = value)) +
geom_histogram(aes(fill = name),
show.legend = FALSE) + # 凡例を表示しない
facet_wrap(vars(as_factor(name))) # 登場順に横に並べる 
一様分布においては，最大最小の平均が良い推定量であることがわかる．
分位点を使ってもう少し詳しくみてみる．
mc_data |> summary() # 四分位点を表示 (簡便な方法) xbar med mid
Min. : 9.445 Min. : 4.318 Min. :11.25
1st Qu.:21.861 1st Qu.:20.093 1st Qu.:23.31
Median :24.916 Median :25.051 Median :24.99
Mean :24.969 Mean :25.024 Mean :24.98
3rd Qu.:28.061 3rd Qu.:29.819 3rd Qu.:26.65
Max. :40.656 Max. :44.534 Max. :37.69 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 4.32 20.1 25.1 29.8 44.5
2 mid 11.3 23.3 25.0 26.6 37.7
3 xbar 9.44 21.9 24.9 28.1 40.7mc_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 4.32 0
2 med 20.1 0.25
3 med 25.1 0.5
4 med 29.8 0.75
5 med 44.5 1
6 mid 11.3 0
7 mid 23.3 0.25
8 mid 25.0 0.5
9 mid 26.6 0.75
10 mid 37.7 1
11 xbar 9.44 0
12 xbar 21.9 0.25
13 xbar 24.9 0.5
14 xbar 28.1 0.75
15 xbar 40.7 1 それぞれのヒストグラムを描くこともできる．
for(name in c("xbar","med","mid")) {
gg <-
mc_data |>
ggplot(aes(x = !!sym(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)
}
違う分布で試してみる．
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(vars(as_factor(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) の項目について以下の問に答えよ．
- 全データを用いてヒストグラム(密度)を作成しなさい．
- ガンマ分布でモデル化して最尤推定を行いなさい．
- 推定した結果をヒストグラムに描き加えて比較しなさい．
- 余力のある者は，自身で収集したデータを用いてモデル化と最尤推定を試みよ．
ヒント
ガンマ分布の最尤推定量は以下のようにして作成できる
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)) # 推定値のみ返す
}解答欄
解答例
データを読み込む．
tw_data <- read_csv("data/tokyo_weather.csv")風速データのヒストグラムを描く．
tw_data |>
ggplot(aes(x = wind)) +
geom_histogram(aes(y = after_stat(density)), # 密度表示
bins = 25, # ビンの数の指定
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))) # 対数密度
## 最尤推定(負の尤度の最小化)の実行
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 <-
tw_data |>
with(mle.gamma(wind)) # データフレームの列名を関数に渡す
print(theta) nu alpha
12.624666 4.608925 得られた結果を前出のヒストグラム上に重ね描きする．
last_plot() +
geom_function(fun = dgamma,
args = list(shape = theta[1],
rate = theta[2]),
colour = "orange",
linewidth = 1.2)
シミュレーションにより推定量の良さの検証(一致性や不偏性の確認)を行うと以下のようになる．
set.seed(5678) # 乱数のシード値の指定 (変更して再現性を確認してみよ)
nu <- 5; alpha <- 2 # 真のパラメタ
grand_truth <- # 真のパラメタのデータフレーム
tibble(nu,alpha) |>
pivot_longer(everything())
mc <- 1000 # 実験回数 (計算が重いので少なめにしている)
#' Monte-Carlo実験
mc_data <-
tibble(n = rep(c(10,50,100), each = mc)) |> # 実験設定(データ数)のデータフレーム
rowwise() |> # 各行ごとに次の処理を実行
mutate(est = list(as_tibble_row( # 実験結果をtibble型式に変換して保存
mle.gamma(rgamma(n,nu,alpha))))) |>
unnest(est) # リストを展開して通常のデータフレームに変換
mc_data |>
mutate(n = as_factor(n)) |> # 実験設定を因子化して整理
pivot_longer(!n) |>
ggplot(aes(x = value, fill = name)) + # nu/alphaごとに密度を推定
geom_density(alpha = 0.6,
show.legend = FALSE) +
geom_vline(data = grand_truth, # 真値の表示
aes(xintercept = value),
colour = "tomato") +
facet_grid(rows = vars(n), cols = vars(name), # nu/alpha およびデータ数ごとにグラフを表示
scales = "free_x") # 座標軸は各列で適宜合わせる
個別に描く場合は，例えば以下のようにすればよい．
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回繰り返した際，真の値(全データによる平均値)が信頼区間に何回含まれるか確認しなさい．
- 余力のある者は，自身で収集したデータで区間推定を試みよ．
解答欄
解答例
全データによる平均値の計算を行う．
mu <-
tw_data |>
pull(solar) |>
mean() # 真値(全データによる平均値)
print(mu)[1] 14.28529これを真値として信頼区間の精度の検証を行う．
ランダムに選んだデータ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
14.30956 17.84124 これを複数回行い，信頼区間の正答率を評価する．
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))) # 信頼区間
}
#' 信頼区間のMonte-Carlo実験
mc_data <-
replicate(mc, mc_trial(n)) |>
t() |> as_tibble() |>
mutate(flag = L < mu & mu < U) # 真値が信頼区間に含まれるか
print(mc_data) # 実験結果を見る# A tibble: 100 × 3
L U flag
<dbl> <dbl> <lgl>
1 11.9 15.5 TRUE
2 12.9 16.7 TRUE
3 12.6 16.5 TRUE
4 11.6 14.6 TRUE
5 11.5 15.0 TRUE
6 13.3 17.1 TRUE
7 13.7 17.3 TRUE
8 11.7 15.3 TRUE
9 13.1 16.2 TRUE
10 13.0 16.0 TRUE
# ℹ 90 more rowsmc_data |> count(flag) # 頻度を見る# A tibble: 2 × 2
flag n
<lgl> <int>
1 FALSE 6
2 TRUE 94信頼区間について更に多数で評価してみる．
mc <- 2000
mc_data <- # 信頼区間のMonte-Carlo実験
replicate(mc, mc_trial(n)) |> t() |> as_tibble() |>
mutate(flag = L < mu & mu < U)
mc_data |>
count(flag) |> # 正否の頻度を見る
mutate(prob = n/sum(n)) # 確率を見る# A tibble: 2 × 3
flag n prob
<lgl> <int> <dbl>
1 FALSE 153 0.0765
2 TRUE 1847 0.924 日射量の平均の推定における観測データを全て使った推定値(真値に相当)と信頼区間の関係をグラフを描いて確認してみる．
k <- 40
mc_data |>
slice_sample(n = k) |> # k個ランダムに選んで描く
rowid_to_column(var = "index") |> # 抜き出したデータに番号を振る
mutate(flag = L < mu & mu < U) |> # 区間に真値が含まれるか判別
ggplot(aes(x = index)) +
geom_errorbar(aes(ymin = L, # 下側信頼限界
ymax = U, # 上側信頼限界
colour = flag), # 判別結果で色分け
width = 0.2) +
geom_hline(yintercept = mu, # 観測データを全て使った推定値
colour = "orange",
linewidth = 1.1)