統計データ解析I

第9講練習問題解答例

Published

June 21, 2025

準備

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

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))}

平均・分散・標準偏差の計算

問題

東京の気候データ (tokyo_weather.csv) の中の

気温 (temp) ，
日射量 (solar) ，
風速 (wind)

の項目について以下の問に答えよ．

全てのデータを用いて各項目の平均・分散・標準偏差を求めよ．(データ数365)
毎月5日のデータのみを用いて各項目の平均・分散・標準偏差を求めよ．(データ数12)
5の付く日(各月の5,15,25)のデータを用いて各項目の平均・分散・標準偏差を求めよ．(データ数36)
ランダムに選んだ36日分のデータで各項目の平均・分散・標準偏差を求めたとき，推定量のばらつきを確認せよ．

Note

データの読み込みは以下のようにすればよい．

tw_data <- read_csv("data/tokyo_weather.csv") # 読み込み方の例

解答例

データの読み込みを行う．

tw_data <- read_csv("data/tokyo_weather.csv")

データフレームの特定の列を取り出し，統計量を計算する愚直な方法は，例えば気温の平均値を計算する場合は以下のようになる．

tw_data |>
    pull(temp) |> # temp列をベクトルとして取り出す
    mean()        # 平均を計算する

[1] 17.64809

列の抽出と計算を同時に行うには関数 dplyr::summarise() を用いれば良い．このとき出力はtibble型のデータフレームになる．

tw_data |>
    summarize(temp_mean = mean(temp)) # temp_mean列が作成される

# A tibble: 1 × 1
  temp_mean
      <dbl>
1      17.6

複数の列に渡って計算する場合には，関数 dplyr::across() を利用する．

tw_data |>
    summarise(across(c(temp,solar,wind), mean))

# A tibble: 1 × 3
   temp solar  wind
  <dbl> <dbl> <dbl>
1  17.6  13.9  2.72

tw_data |>
    summarise(across(c(temp,solar,wind), var))

# A tibble: 1 × 3
   temp solar  wind
  <dbl> <dbl> <dbl>
1  69.9  49.3 0.788

tw_data |>
    summarise(across(c(temp,solar,wind), sd))

# A tibble: 1 × 3
   temp solar  wind
  <dbl> <dbl> <dbl>
1  8.36  7.02 0.888

統計量をまとめて計算するには以下のようにすれば良い．集計されたデータフレームには “データの列名”_“関数のリスト名”の列が組合せの数だけ作成される．

(tw_summary <- # 後で参照するために保存しておく
     tw_data |> 
     summarise(across(c(temp,solar,wind),
                      list(mean = mean, var = var, sd = sd))))

# A tibble: 1 × 9
  temp_mean temp_var temp_sd solar_mean solar_var solar_sd wind_mean wind_var
      <dbl>    <dbl>   <dbl>      <dbl>     <dbl>    <dbl>     <dbl>    <dbl>
1      17.6     69.9    8.36       13.9      49.3     7.02      2.72    0.788
# ℹ 1 more variable: wind_sd <dbl>

データの部分集合，例えば毎月5日のデータによる計算は関数 dplyr::filter() を用いて以下のようにすればよい．

tw_data |>
    filter(day == 5) |>
    summarise(across(c(temp,solar,wind),
                     list(mean = mean, var = var, sd = sd)))

# A tibble: 1 × 9
  temp_mean temp_var temp_sd solar_mean solar_var solar_sd wind_mean wind_var
      <dbl>    <dbl>   <dbl>      <dbl>     <dbl>    <dbl>     <dbl>    <dbl>
1      17.3     87.0    9.33       13.8      116.     10.8      2.48    0.312
# ℹ 1 more variable: wind_sd <dbl>

同様に5の付く日のデータによる計算は以下のようになる．

tw_data |>
    filter(day %in% c(5,15,25)) |>
    summarise(across(c(temp,solar,wind),
                     list(mean = mean, var = var, sd = sd)))

# A tibble: 1 × 9
  temp_mean temp_var temp_sd solar_mean solar_var solar_sd wind_mean wind_var
      <dbl>    <dbl>   <dbl>      <dbl>     <dbl>    <dbl>     <dbl>    <dbl>
1      17.7     76.0    8.72       14.5      76.0     8.72      2.67    0.535
# ℹ 1 more variable: wind_sd <dbl>

ランダムに選択した36日で推定した場合のばらつきを調べる．行をランダムに選択するには関数 dplyr::slice_sample() が利用できる．

mc <- 5000 # 実験回数を指定
n <- 36    # ランダムに選択する日数を指定
#' 気温の標本平均による例
my_trial <- function(){ 
    tw_data |>
        slice_sample(n = n) |>   # ランダムにn行抽出
        summarise(mean(temp)) |> # tibble形式
        pull()                   # ベクトル形式(単なる数値)に変換
}
xbars <- replicate(mc, my_trial())
tibble(平均気温の推定値 = xbars) |>
    ggplot(aes(x = 平均気温の推定値)) +
    geom_histogram(aes(y = after_stat(density)), 
                   bins = 30,                    
                   fill ="lightblue",          
                   colour = "blue") + 
    geom_vline(xintercept = tw_summary[["temp_mean"]], # 全体の平均
               colour = "red")

Note

上記の例では保存しておいた平均の値を参照しているが， geom_vline の中で計算し直すこともできる．

tibble(平均気温の推定値 = xbars) |>
    ggplot(aes(x = 平均気温の推定値)) +
    geom_histogram(aes(y = after_stat(density)), 
                   bins = 30,                    
                   fill ="lightblue",          
                   colour = "blue") + 
    geom_vline(xintercept = tw_data |> 
                   summarise(mean(temp)) |>
                   pull(), # tibble形式を数値に変換
               colour = "red")

Note

例えば以下のようにすれば，項目および統計量を総当たりで調べることができる．

my_trial <- function(n){ # 日数を指定してまとめて計算するように変更する
    tw_data |>
        slice_sample(n = n) |>
        summarise(across(c(temp,solar,wind),
                         list(mean = mean, var = var, sd = sd))) |>
        unlist() # データフレームではなくベクトルとして返す
}
my_data <- # 実験データをデータフレームに直しておく
    replicate(mc, my_trial(n = n)) |> t() |> as_tibble()
#' 図示する部分を関数として定義しておく
#' 文字列を利用して関数aes()などで列名を指定するには
#'   aes(x = !!sym("文字列"))
#' という非標準評価(NSE)の仕組みを使う必要があるので注意する
my_plot <- function(data,    # データ
                    summary, # データ全体から計算した真値
                    item,    # 項目名 (文字列で指定)
                    func){   # 集計関数名 (文字列で指定)
    name <- sym(paste(item, func, sep = "_")) # シンボルを作成
    gg <- data |>
        ggplot(aes(x = !!name)) + # シンボルを !! によって unquote
        geom_histogram(aes(y = after_stat(density)), 
                       bins = 30,                    
                       fill ="lightblue",          
                       colour = "blue") + 
        geom_vline(xintercept = summary[[name]],
                   colour = "red") +
        labs(title = paste(item, "の", func, "の推定"))
    print(gg) # 関数内では明示的に描画を指示する
}
#' 全ての組み合わせは for 文で実行可能
for(item in c("temp","solar","wind")){ # 項目を指定
    for(func in c("mean","var","sd")){ # 関数を指定
        my_plot(my_data, tw_summary, item, func)
    }
}

いずれの変数においても標準偏差の推定量は若干の偏りがあることが確認できる．

風速の分散の推定量の分布が他と異なり正規分布に近くないので，サンプル数を3倍に増やしてみる．

my_plot(replicate(mc, my_trial(n = 108)) |> t() |> as_tibble(), 
        tw_summary, "wind", "var")

推定量の分散が小さくなるとともに形状が正規分布に近づいたことが確認できる．

歪度と超過尖度の計算

問題

東京の気候データ (tokyo_weather.csv) の中の

気温 (temp) ，
日射量 (solar) ，
風速 (wind)

の項目について以下の問に答えよ．

全てのデータを用いて各項目の歪度と超過尖度を求めよ．(データ数365)
5のつく日のデータのみを用いて各項目の歪度と超過尖度を求めよ．(データ数36)
それぞれの値から正規分布から逸脱していると思われる項目はいずれか考察せよ．
各データのヒストグラムを描き，データから計算される平均と分散を持つ正規分布と比較せよ．

解答例

歪度と超過尖度を計算するためのパッケージを読み込む．

library("e1071")

全データによる計算は以下のようにすれば良い．

tw_data |> 
    summarise(across(c(temp,solar,wind),
                     list(skew = skewness, kurt = kurtosis)))

# A tibble: 1 × 6
  temp_skew temp_kurt solar_skew solar_kurt wind_skew wind_kurt
      <dbl>     <dbl>      <dbl>      <dbl>     <dbl>     <dbl>
1   -0.0109     -1.32      0.246     -0.796      1.41      2.86

Note

パッケージの一部の関数しか利用しない場合は，関数 library() でパッケージを読み込まなくても，パッケージ名とともに関数を指定することで利用することができる．ただし，関数内で他の関数に依存する場合は注意する必要がある．

tw_data |> 
    summarise(across(c(temp,solar,wind),
                     list(skew = e1071::skewness,
                          kurt = e1071::kurtosis)))

# A tibble: 1 × 6
  temp_skew temp_kurt solar_skew solar_kurt wind_skew wind_kurt
      <dbl>     <dbl>      <dbl>      <dbl>     <dbl>     <dbl>
1   -0.0109     -1.32      0.246     -0.796      1.41      2.86

前問と同様に，5の付く日のデータによる計算は以下のようになる．

tw_data |>
    filter(day %in% c(5,15,25)) |>
    summarise(across(c(temp,solar,wind),
                     list(skew = skewness, kurt = kurtosis)))

# A tibble: 1 × 6
  temp_skew temp_kurt solar_skew solar_kurt wind_skew wind_kurt
      <dbl>     <dbl>      <dbl>      <dbl>     <dbl>     <dbl>
1    -0.104     -1.33     0.0314      -1.40     0.930     0.789

Note

歪度や尖度は3次・4次の計算なので，少数のデータでは計算結果が不安定になりやすい．このためデータ数については注意する必要がある．

各データのヒストグラムと正規分布を比較すると以下のようになる．

tw_summary <- # 前の練習問題の結果を利用 
    tw_data |> 
    summarise(across(c(temp,solar,wind),
                     list(mean = mean, var = var, sd = sd)))
for(item in c("temp","solar","wind")){ # 項目を指定
    item <- sym(item)
    gg <- tw_data |>
        ggplot(aes(x = !!item)) + 
        geom_histogram(aes(y = after_stat(density)), 
                       bins = 30,                    
                       fill ="lightblue",          
                       colour = "blue") + 
        geom_function(fun = dnorm, # 標本平均と標本標準偏差を利用する
                      args = list(mean = tw_summary[[paste0(item,"_mean")]],
                                  sd = tw_summary[[paste0(item,"_sd")]]),
                      colour = "red") +
        labs(title = paste(item, "と正規分布の比較"))
    print(gg)
}

共分散と相関の計算

問題

東京の気候データ (tokyo_weather.csv) の中の

気温 (temp)
降水量 (rain)
日射量 (solar)
風速 (wind)
気圧 (press)
湿度 (humid)

の項目(いずれも数値データ)について以下の問に答えよ．

それぞれの項目間の共分散，および相関を求めよ．
相関の高い項目の組(絶対値が大きい)，および相関の低い項目の組(0に近い)を求めよ．
その項目同士の散布図を描け．

解答例

共分散および相関行列の計算を行う．

tw_cov <-
    tw_data |>
    select(temp,rain,solar,wind,press,humid) |>
    cov()
tw_cor <-
    tw_data |>
    select(temp,rain,solar,wind,press,humid) |>
    cor()

特徴的な2変数間の関係を調べてみる．

tw_cor == min(tw_cor)                    # 負の最大相関 (temp,press)

       temp  rain solar  wind press humid
temp  FALSE FALSE FALSE FALSE  TRUE FALSE
rain  FALSE FALSE FALSE FALSE FALSE FALSE
solar FALSE FALSE FALSE FALSE FALSE FALSE
wind  FALSE FALSE FALSE FALSE FALSE FALSE
press  TRUE FALSE FALSE FALSE FALSE FALSE
humid FALSE FALSE FALSE FALSE FALSE FALSE

tw_cor == max(tw_cor-diag(diag(tw_cor))) # 対角を除く最大相関 (temp,humid)

       temp  rain solar  wind press humid
temp  FALSE FALSE FALSE FALSE FALSE  TRUE
rain  FALSE FALSE FALSE FALSE FALSE FALSE
solar FALSE FALSE FALSE FALSE FALSE FALSE
wind  FALSE FALSE FALSE FALSE FALSE FALSE
press FALSE FALSE FALSE FALSE FALSE FALSE
humid  TRUE FALSE FALSE FALSE FALSE FALSE

abs(tw_cor) == min(abs(tw_cor))          # 最小相関(0に近い) (rain,wind)

       temp  rain solar  wind press humid
temp  FALSE FALSE FALSE FALSE FALSE FALSE
rain  FALSE FALSE FALSE  TRUE FALSE FALSE
solar FALSE FALSE FALSE FALSE FALSE FALSE
wind  FALSE  TRUE FALSE FALSE FALSE FALSE
press FALSE FALSE FALSE FALSE FALSE FALSE
humid FALSE FALSE FALSE FALSE FALSE FALSE

散布図を描画する．

library(GGally) # 散布図行列を描くためのパッケージ
tw_data |> # 対象データを全て表示してみる
    select(temp,rain,solar,wind,press,humid) |>
    ggpairs()

tw_data |> # columnsを指定すれば関数selectを使わなくてもよい
    ggpairs(columns = c("temp","rain","solar","wind","press","humid"))

tw_data |> # 月ごとに色を変えることもできる
    ggpairs(columns = c("temp","rain","solar","wind","press","humid"),
            lower = list(mapping = aes(colour = as_factor(month))))

tw_data |> # 負の最大相関
    ggpairs(columns = c("temp","press"))

tw_data |> # 最大相関
    ggpairs(columns = c("temp","humid"))

tw_data |> # 最小相関
    ggpairs(columns = c("rain","wind"))

Note

数値項目を抽出するのであれば以下のように簡潔に書ける．

tw_data |>
    select(where(is.numeric))

# A tibble: 366 × 11
    year month   day  temp  rain solar  snow  wind press humid cloud
   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1  2024     1     1   9.2     0 11.7      0   4   1016.    48   2  
 2  2024     1     2   6.1     0  5.43     0   1.8 1017.    65   6.8
 3  2024     1     3   7.4     0  5.91     0   1.5 1010.    72   7.5
 4  2024     1     4   9.7     0 10.3      0   3.1 1010.    58   4.8
 5  2024     1     5   7.4     0 11.4      0   1.7 1016     57   0  
 6  2024     1     6   9.1     0 11.4      0   1.4 1009.    61   0  
 7  2024     1     7   8.3     0  7.85     0   2.1 1005.    53   6.3
 8  2024     1     8   5.4     0 12.3      0   3.4 1015.    43   0  
 9  2024     1     9   5.4     0 11.9      0   1.8 1014.    49   1.8
10  2024     1    10   7.2     0  9.35     0   1.8 1010.    52   5  
# ℹ 356 more rows

ただし不要な項目も含まれるので，別途以下のような条件を追加するなど工夫は必要となる．

tw_data |>
    select(where(is.numeric) & !year:day & !c(snow,cloud))

# A tibble: 366 × 6
    temp  rain solar  wind press humid
   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1   9.2     0 11.7    4   1016.    48
 2   6.1     0  5.43   1.8 1017.    65
 3   7.4     0  5.91   1.5 1010.    72
 4   9.7     0 10.3    3.1 1010.    58
 5   7.4     0 11.4    1.7 1016     57
 6   9.1     0 11.4    1.4 1009.    61
 7   8.3     0  7.85   2.1 1005.    53
 8   5.4     0 12.3    3.4 1015.    43
 9   5.4     0 11.9    1.8 1014.    49
10   7.2     0  9.35   1.8 1010.    52
# ℹ 356 more rows

分位点と最頻値の計算

問題

東京の気候データ (tokyo_weather.csv) の中の

気温 (temp; 数値データ)
最多風向 (wdir; ラベルデータ)

を用いて以下の問に答えよ．

全てのデータを用いて気温の四分位点を求めよ．(データ数365)
5の付く日(各月の5,15,25)の気温の四分位点を求めよ．(データ数36)
ランダムに選んだ36日分のデータで気温の四分位点がどのくらいばらつくか確認せよ．
風向の最頻値を求めよ．

解答例

気温の分位点を求める．全データによる計算は以下のようになる．

(tw_temp_summary <-
     tw_data |>
     pull(temp) |> # 1列のみ抽出してベクトルにする
     summary())

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.800   9.425  18.300  17.648  24.600  32.300

データの部分集合による計算(例えば5の付く日のデータ)は以下のようにすればよい．

tw_data |>
    filter(day %in% c(5,15,25)) |>
    pull(temp) |>
    summary()

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   3.00   10.95   19.45   17.74   24.65   30.30

ランダムに選択した36日で推定した場合のばらつきを調べる．

mc <- 5000 # 実験回数を指定
my_trial <- function(){
    tw_data |>
        slice(sample(nrow(tw_data),36)) |> # 重複なしで36行選ぶ
        pull(temp) |>
        summary()
}
my_data <-
    replicate(mc, my_trial()) |> t() |> as_tibble()
#' ヒストグラムの表示
for(name in names(my_data)){
    name <- sym(name)
    gg <- my_data |>
        ggplot(aes(x = !!name)) +
        geom_histogram(aes(y = after_stat(density)),
                       bins = 40,
                       fill = "lightblue",
                       colour = "blue") +
        geom_vline(xintercept = tw_temp_summary[[name]],
                   colour = "red") +
        labs(x = "気温",
             title = paste("気温の", name, "の推定のばらつき"))
    print(gg)
}

以下のようにすれば定義域とビンを揃えて表示することができる．

breaks <- # 全データを用いて適切なビンの計算して固定する
    tw_data |> pull(temp) |> pretty(n = 50)
for(name in names(my_data)){
    name <- sym(name)
    gg <- my_data |>
        ggplot(aes(x = !!name)) +
        geom_histogram(aes(y = after_stat(density)),
                       breaks = breaks, # ビンを固定
                       fill = "lightblue",
                       colour = "blue") +
        geom_vline(xintercept = tw_temp_summary[[name]],
                   colour = "red") +
        ylim(0,0.5) + # y軸の表示を適当に固定する
        labs(x = "気温",
             title = paste("気温の", name, "の推定のばらつき"))
    print(gg)
}

表示には密度推定を用いても良い．

for(name in names(my_data)){
    name <- sym(name)
    gg <- my_data |>
        ggplot(aes(x = !!name)) +
        geom_density(fill = "lightblue",
                     colour = "blue") +
        geom_vline(xintercept = tw_temp_summary[[name]],
                   colour = "red") +
        lims(x = range(tw_data |> pull(temp)), y = c(0,0.5)) +
        labs(x = "気温",
             title = paste("気温の", name, "の推定のばらつき"))
    print(gg)
}

最多風向の最頻値を調べる．

(my_table <-
     tw_data |>
     pull(wdir) |>
     table()) # 頻度表


  E ENE ESE   N  NE NNE NNW  NW   S  SE SSE SSW  SW WNW 
  2  13   1  11  18  22  57  73  56   7  79  15   1  11

max(my_table) # 最大値

[1] 79

which.max(my_table) # 最頻値の表の位置

SSE 
 11

names(which.max(my_table)) # 最頻値の項目名

[1] "SSE"

Note

例えば以下のようにすれば16方位の頻度をグラフ化することができる．

cardinal_directions <- c("N","NNE","NE","ENE",
                         "E","ESE","SE","SSE",
                         "S","SSW","SW","WSW",
                         "W","WNW","NW","NNW")
tw_data |>
    mutate(wdir = factor(wdir, levels = cardinal_directions)) |>
    group_by(wdir, .drop = FALSE) |> # 0も含めて16方位全て計算する
    summarize(count = n()) |> # 関数 table() だと頻度 0 の項目はないことに注意
    ggplot(aes(x = wdir, y = count)) +
    geom_bar(stat = "identity", fill = alpha("blue", 0.4)) +
    coord_polar(start = 0) # 極座標表示にする

中心を空洞にするには負の値で下駄を履かせればよい．

tw_data |>
    mutate(wdir = factor(wdir, levels = cardinal_directions)) |>
    group_by(wdir, .drop = FALSE) |> # 0も含めて16方位全て計算する
    summarize(count = n()) |>
    ggplot(aes(x = wdir, y = count)) +
    geom_bar(stat = "identity", fill = alpha("blue", 0.4)) +
    ylim(-40,80) + # 負の側を表示範囲に加える
    coord_polar(start = -pi/16) + # 北を真上に配置
    labs(x = NULL, y = NULL, title = "風向きの頻度") # 不要なラベルを削除