---
title: "統計データ解析I"
subtitle: "第8講 練習問題 解答例"
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_{1},\dotsc,X_n$ を成功確率 $p$ の Bernoulli 分布に従う独立同分布な確率変数列とする． このとき $\sum_{i=1}^nX_i$ の分布は，試行回数 $n$ , 成功確率 $p$ の二項分布に従う．これをグラフに描画して確認しなさい．

### 解答例

試行を模擬する関数を定義する．

```{r}
n <- 16 
p <- 0.6
my_random <- function(){ # Bernolli分布をm個生成して合計
    sum(rbinom(n, size=1, prob=p))}
```

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

```{r}
mc <- 10000 # 実験回数を指定
my_data <- replicate(mc, my_random())
```

出現値ごとの度数分布表，および正規化した確率表を表示する．

```{r}
my_data |> table()                 # base::table による頻度表の作成
my_data |> table() |> prop.table() # 確率値に変換
```

`tibble`形式のデータフレームを作成するには関数 `dplyr::count()` を利用して例えば以下のようにすればよい．

```{r}
my_data |>
    as_tibble_col() |> # ベクトルから1列のtibbleを作成．列名の既定値は value
    count(value)       # value 列の頻度を集計する
```

棒グラフを利用して頻度表を視覚化する．

```{r}
my_data |> as_tibble_col() |>
    count(value) |>
    ggplot(aes(x = value, y = n)) +
    geom_bar(stat = "identity") # 最も単純な表示
```

理論値(確率)は関数 `dbinom()` で計算できるので，組み合わせて横並びに表示する．

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

::: callout-note
上記の例では `dplyr::count` を利用する場合を示したが， 視覚化のためのデータフレームはいろいろな方法で作成できる．

`base::table()` および `balse::prop.table()` を利用する場合．

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

`package::janitor` を利用する場合．

```{r}
library(janitor)
my_data |> as_tibble_col() |> 
    tabyl(value) |>                                # 確率
    set_names(c("値","頻度","観測値")) |>            # 列名を変更
    mutate(理論値 = dbinom(値, size = n, prob = p)) # 理論値を追加
```
:::

## Poisson 分布

### 問題

Poisson 分布に関して以下を考察せよ．

-   平均と分散の理論計算を確認しなさい．

-   $X,Y$ を独立な2つの確率変数とし，それぞれ強度 $\lambda_{1},\lambda_{2}$ の Poisson 分布に従うとする． このとき, 和 $X+Y$ の分布は強度 $\lambda_{1}+\lambda_{2}$ の Poisson 分布に従う． これをグラフに描画して確認しなさい．

### 解答例

試行を模擬する関数を定義する．

```{r}
lambda1 <- 5
lambda2 <- 12
my_random <- function(){ # 2つの Poisson 分布の和
    rpois(1, lambda=lambda1)+rpois(1, lambda=lambda2)}
```

確率シミュレーションを行ない，視覚化する．

```{r}
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}$ は独立かつともに標準正規分布に従う． これをグラフに描画して確認しなさい．

    ::: callout-note
    この変換を Box-Muller 変換と呼ぶ．
    :::

### 解答例

試行を模擬する関数を作成する．

```{r}
my_random <- function(){ # 一方の分布を確認する
    u1 <- runif(1)
    u2 <- runif(1)
    return(sqrt(-2*log(u1))*cos(2*pi*u2))}
```

確率シミュレーションを行い，ヒストグラムを用いて視覚化する．

```{r}
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()` を用いて重ね描きする．

```{r}
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つの確率変数の関係を調べてみる． それぞれを取り出せるように関数を作成する．

```{r}
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つの変数の散布図を描く．

```{r}
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,]` などとすれば良い．

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

::: callout-note
散布図を描くためのデータフレームは以下のようにして作ることもできる． 高次元の場合には関数を利用した方が簡単に記述できる．

```{r}
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$ 分布に従うことを確認しなさい．

### 解答例

試行を模擬する関数を定義する．

```{r}
k <- 8 # 自由度を設定
my_random <- function(){ 
    sum(rnorm(k)^2)} # k個の標準正規乱数の二乗和
```

確率シミュレーションを行い，視覚化する．

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

::: callout-note
関数 bquote() は関数 expression() と同様な働きをするが，自由度の k を .(k) で評価して取り込むことができる． 以下と表示を比較してみよう．

```{r}
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$ 分布に従うことを確認しなさい．

### 解答例

試行を模擬する関数を作成する．

```{r}
k <- 8 # 自由度を設定
my_random <- function(){ 
    sum(rnorm(k)^2)} # k個の標準正規乱数の二乗和
```

確率シミュレーションを行い，視覚化する．

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

::: callout-note
関数 bquote() は関数 expression() と同様な働きをするが， 自由度の k を .(k) で評価して取り込むことができる． 以下と比較してみよう．

```{r}
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$ 分布に従うことを確認しなさい．

### 解答例

試行を模擬する関数を作成する．

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

確率シミューレションを行い，視覚化する．

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

::: callout-note
グラフの一部に着目したい場合は xlim/ylim で表示を調整することができる．

```{r}
last_plot() +               # 直前のプロット
    xlim(0,6) + ylim(0,0.9) # 表示領域を指定する(表示しない部分について警告が出る)
```
:::