---
title: "統計データ解析I"
subtitle: "第10講 練習問題 解答例"
date: "`r Sys.time()`"
format:
    html: 
      toc: true
      html-math-method: katex
      self-contained: true
      grid: 
        margin-width: 350px
execute: 
  echo: true
  warning: false
reference-location: margin
citation-location: margin
tbl-cap-location: margin
fig-cap-location: margin
editor: visual
editor_options: 
  chunk_output_type: console
---

## 準備

以下で利用する共通パッケージを読み込む．

```{r}
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))}
```

## 一様分布の平均の推定

### 問題

$X$ を一様乱数に従う確率変数とし，平均値の推定量として以下を考える． それぞれの推定量の分散を比較しなさい．

-   標本平均 (mean)
-   中央値 (median)
-   最大値と最小値の平均 ((max+min)/2)

::: callout-note
#### ヒント

以下のような関数を作り，確率シミュレーションを行えばよい．

```{r}
#| eval: false
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つまとめて計算する関数
```
:::

### 解答例

平均値の推定を行う関数(標本平均，中央値，最大最小の平均)を定義する．

```{r}
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)) # 最大最小の平均
}
```

実験の設定を宣言する．

```{r}
n <- 10         # 観測データのサンプル数
mc <- 10000     # 実験回数
a <- 0; b <- 50 # 一様乱数の区間
```

確率シミュレーションを行う．

```{r}
mc_data <- 
    replicate(mc, estimate_means(n, min = a, max = b)) |> # 3 x mc 型行列
    t() |>      # 転置して mc x 3 型行列に変換
    as_tibble() # データフレーム化
mc_data # 実験結果の一部を確認
```

それぞれの推定量の平均と分散を計算し，結果を集計する．

```{r}
mc_data |> 
    summarize(across(everything(),
                     list(mean = mean, var = var)))
```

別々のヒストグラムにして比較してみる．

```{r}
mc_data |>
    pivot_longer(everything()) |>
    ggplot(aes(x = value)) +
    geom_histogram(aes(fill = name),
                   show.legend = FALSE) + # 凡例を表示しない
    facet_wrap(~name)                     # 横に並べる 
```

一様分布においては，最大最小の平均が良い推定量であることがわかる．

::: callout-note
分位点を使ってもう少し詳しくみてみる．

```{r}
mc_data |> summary() # 四分位点を表示 (簡便な方法)
mc_data |> # データフレームとして取得する方法 (以下は一例)
    pivot_longer(everything()) |> # long format に変更
    group_by(name) |> # 推定量ごとに集計
    summarize(as_tibble_row(quantile(value))) # 計算結果を行に並べる
mc_data |> # 別形式のデータフレームの作成例
    pivot_longer(everything()) |> # long format に変更
    group_by(name) |> # 推定量ごとに集計
    reframe(qs=quantile(value), probs=seq(0,1,.25)) 
```

それぞれのヒストグラムを描くこともできる．

```{r}
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)
}
```
:::

::: callout-note
違う分布で試してみる．

```{r}
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) 
```

この分布では中央値が良い推定量となることがわかる．
:::

::: callout-tip
Rで用いることのできる色の名前は関数 colors()/colours() で確認することができる. 例えば以下のようにすると視覚的に確認することができる.

```{r}
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`) の項目について以下の問に答えよ．

-   全データを用いてヒストグラム(密度)を作成しなさい．
-   ガンマ分布でモデル化して最尤推定を行いなさい．
-   推定した結果をヒストグラムに描き加えて比較しなさい．

::: callout-note
ガンマ分布の最尤推定量は以下のようにして作成できる

```{r}
#| eval: false
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)) # 推定値のみ返す
}
```
:::

::: callout-note
余力のある者は，自身で収集したデータを用いてモデル化と最尤推定を試みよ．
:::

### 解答例

データを読み込む．

```{r}
tw_data <- read_csv("data/tokyo_weather.csv")
```

風速データのヒストグラムを描く．

```{r}
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 = "風速のヒストグラム")
```

ガンマ分布の最尤推定量を求めるための関数を作成する．

```{r}
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)) # 推定値のみ返す
    }
}
```

最尤推定を行う．

```{r}
(theta <- with(tw_data, mle.gamma(wind))) 
```

得られた結果を前出のヒストグラム上に重ね描きする．

```{r}
last_plot() +
    geom_function(fun = dgamma,
                  args = list(shape = theta[1],
                              rate = theta[2]),
                  colour = "orange",
                  linewidth = 1.2)
```

::: callout-note
シミュレーションにより推定量の良さの検証(一致性や不偏性の確認)を行うと以下のようになる．

```{r}
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回繰り返した際，真の値(全データによる平均値)が信頼区間に何回含まれるか確認しなさい．

::: callout-note
余力のある者は，自身で収集したデータで区間推定を試みよ．
:::

### 解答例

全データによる平均値の計算を行う．

```{r}
(mu <- tw_data |> pull(solar) |> mean()) # 真値
```

これを真値として信頼区間の精度の検証を行う．

ランダムに選んだデータ50点を用いて90%信頼区間を構成する．

```{r}
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))) # 信頼区間
```

これを複数回行い，信頼区間の正答率を評価する．

```{r}
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)             # 頻度を見る
```

::: callout-note
信頼区間について更に多数で評価してみる．

```{r}
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)) # 確率を見る
```
:::

::: callout-note
日射量の平均の推定における観測データを全て使った推定値(真値に相当)と信頼区間の関係をグラフを描いて確認してみる．

```{r}
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)
```
:::