統計データ解析

第2講 サンプルコード

準備

データの操作と視覚化のために tidyverse パッケージを読み込む．

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Quartoファイルの中で使われているパッケージがインストールされていない場合は， エディタ上部にインストールを促す警告が出るので，その指示に従ってインストールすればよい． 手動でインストールする場合は以下のいずれかを実行する．

• Package タブからインストール
• コンソール上で次のコマンドを実行 install.packages("パッケージ名")

tidyverse は第1講のサンプルで使っているので通常はインストールされているが， R言語のバージョンがアップグレードされた場合などに対応が必要となる．

library(tidyverse) を実行すると読み込まれたパッケージが表示される． 同名の関数が存在する場合には Conflicts として衝突する関数名が表示される． 衝突する場合はパッケージ名を明示的に付ける必要がある． 例えば filter() には以下の2つがある．

• dplyer::filter() データフレームの抽出のための関数
• stats::filter() 時系列処理のための線形フィルタ関数

Tip

名前の衝突による不具合を避けたい場合は conflicted パッケージの利用を推奨する．

library(conflicted)

データフレームの作成

問題

次の表に対応するデータフレームを作成しなさい．

name	math	phys	chem	bio
Alice	90	25	65	70
Bob	80	50	100	50
Carol	70	75	70	30
Dave	60	100	40	80
Eve	50	80	75	100

解答例

各項目が同じ長さのベクトルなので，これらを並べればよい．

(grade_data <- tibble( # 変数名は自由に決めてよい
     name = c("Alice", "Bob", "Carol", "Dave", "Eve"),
     math = c(90, 80, 70, 60, 50),
     phys = c(25, 50, 75,100, 80),
     chem = c(65,100, 70, 40, 75),
     bio  = c(70, 50, 30, 80,100)))

# A tibble: 5 × 5
  name   math  phys  chem   bio
  <chr> <dbl> <dbl> <dbl> <dbl>
1 Alice    90    25    65    70
2 Bob      80    50   100    50
3 Carol    70    75    70    30
4 Dave     60   100    40    80
5 Eve      50    80    75   100

行や列の名前を操作するには以下のようにする．

names(grade_data)   # 列の名前を表示する

[1] "name" "math" "phys" "chem" "bio" 

print(grade_data)   # データフレームの内容を表示(printは省略できる)

# A tibble: 5 × 5
  name   math  phys  chem   bio
  <chr> <dbl> <dbl> <dbl> <dbl>
1 Alice    90    25    65    70
2 Bob      80    50   100    50
3 Carol    70    75    70    30
4 Dave     60   100    40    80
5 Eve      50    80    75   100

glimpse(grade_data) # データフレームの内容を別の形式で表示

Rows: 5
Columns: 5
$ name <chr> "Alice", "Bob", "Carol", "Dave", "Eve"
$ math <dbl> 90, 80, 70, 60, 50
$ phys <dbl> 25, 50, 75, 100, 80
$ chem <dbl> 65, 100, 70, 40, 75
$ bio  <dbl> 70, 50, 30, 80, 100

データの取り出し方は後ほど詳しく述べるが，例えば以下のように操作することができる．

grade_data[2,3]      # 特定の要素を数値で参照する

# A tibble: 1 × 1
   phys
  <dbl>
1    50

grade_data[2,"phys"] # 列を名前で参照する (上記と同じ結果)

# A tibble: 1 × 1
   phys
  <dbl>
1    50

grade_data[3,]       # 特定の行を表示 (データフレームになる)

# A tibble: 1 × 5
  name   math  phys  chem   bio
  <chr> <dbl> <dbl> <dbl> <dbl>
1 Carol    70    75    70    30

grade_data["bio"]    # 特定の列を表示 (データフレームになる)

# A tibble: 5 × 1
    bio
  <dbl>
1    70
2    50
3    30
4    80
5   100

grade_data[,"bio"]   # 上記と同じ結果

# A tibble: 5 × 1
    bio
  <dbl>
1    70
2    50
3    30
4    80
5   100

grade_data[["bio"]]  # ベクトルとして取り出す (リストとしての処理)

[1]  70  50  30  80 100

grade_data$bio       # 上記と同じ結果

[1]  70  50  30  80 100

Tip

変数や関数の名称については各自で命名規則を決めておくと良い． 例えば以下のサイトが参考になる．

• https://style.tidyverse.org/index.html
• https://google.github.io/styleguide/Rguide.html

ファイルの読み書き

問題

以下の問いに答えなさい．

1. 前の演習で作成したデータフレームを適当なファイルに書き出しなさい．
2. 書き出したファイルから別の変数に読み込みなさい．
3. pcr_case_daily.csv (厚労省からダウンロードしたファイル)を変数 pcr_data に読み込みなさい．

解答例

前の練習問題で作ったデータフレームを作業ディレクトリの data 以下に grade_data.csv として保存する．

write_csv(grade_data, file = "data/grade_data.csv")

保存したファイルは File タブからその中身を含め確認することができる．

grade_copy というオブジェクトにファイルの内容を代入して，その中身を確認する．

(grade_copy <- read_csv(file = "data/grade_data.csv"))

# A tibble: 5 × 5
  name   math  phys  chem   bio
  <chr> <dbl> <dbl> <dbl> <dbl>
1 Alice    90    25    65    70
2 Bob      80    50   100    50
3 Carol    70    75    70    30
4 Dave     60   100    40    80
5 Eve      50    80    75   100

Files タブの操作で読み込みことも可能である．

ダウンロードしたファイルをファイル名 pcr_case_daily.csv として作業ディレクトリの data 以下に保存しておく．

pcr_data <- read_csv("data/pcr_case_daily.csv") 
print(pcr_data) # 中身の一部を表示

# A tibble: 1,175 × 10
   日付    国立感染症研究所 検疫所 地方衛生研究所・保健…¹ 民間検査会社（主に行…²
   <chr>              <dbl>  <dbl>                  <dbl>                  <dbl>
 1 2020/2…              472     75                    398                      0
 2 2020/2…               15     68                    609                      0
 3 2020/2…               20     15                    758                      0
 4 2020/2…              261    188                    902                    132
 5 2020/2…              341    127                    677                      2
 6 2020/2…               53     72                    529                      0
 7 2020/2…               22    103                    460                     17
 8 2020/2…              195     38                    695                      0
 9 2020/2…              267     19                    970                    149
10 2020/2…              237     61                   1101                      0
# ℹ 1,165 more rows
# ℹ abbreviated names: ¹​`地方衛生研究所・保健所`,
#   ²​`民間検査会社（主に行政検査）`
# ℹ 5 more variables: 大学等 <dbl>, 医療機関 <dbl>, 小計 <dbl>,
#   `民間検査会社（主に自費検査）` <dbl>, 計 <dbl>

列名を確認して pcr_colnames に保存しておく．

(pcr_colnames <- names(pcr_data)) # colnames(pcr_data) でも良い

 [1] "日付"                         "国立感染症研究所"            
 [3] "検疫所"                       "地方衛生研究所・保健所"      
 [5] "民間検査会社（主に行政検査）" "大学等"                      
 [7] "医療機関"                     "小計"                        
 [9] "民間検査会社（主に自費検査）" "計"                          

列名を扱い易いように英語略記に変更する．

略称はそれぞれ以下の意味．

• National Institute of Infectious Diseases
• Customs-Immigration-Quarantine
• Health Center
• Administrative Inspection
• University
• Medical Institution
• subtotal
• Self Inspection
• total

names(pcr_data) <- # 
    c("date","niid","ciq","hc","ai","univ","mi","sub","si","total")
pcr_data # 中身を確認(10行だけ表示される)

# A tibble: 1,175 × 10
   date       niid   ciq    hc    ai  univ    mi   sub    si total
   <chr>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1 2020/2/18   472    75   398     0    79    NA  1024    NA  1024
 2 2020/2/19    15    68   609     0     0    NA   692    NA   692
 3 2020/2/20    20    15   758     0     0    NA   793    NA   793
 4 2020/2/21   261   188   902   132   108    NA  1591    NA  1591
 5 2020/2/22   341   127   677     2    19    NA  1166    NA  1166
 6 2020/2/23    53    72   529     0     0    NA   654    NA   654
 7 2020/2/24    22   103   460    17     0    NA   602    NA   602
 8 2020/2/25   195    38   695     0     0    NA   928    NA   928
 9 2020/2/26   267    19   970   149     0    NA  1405    NA  1405
10 2020/2/27   237    61  1101     0     0    NA  1399    NA  1399
# ℹ 1,165 more rows

names(pcr_colnames) <- names(pcr_data) # 和英の列名の対応づけができるようにしておく
pcr_colnames

                          date                           niid 
                        "日付"             "国立感染症研究所" 
                           ciq                             hc 
                      "検疫所"       "地方衛生研究所・保健所" 
                            ai                           univ 
"民間検査会社（主に行政検査）"                       "大学等" 
                            mi                            sub 
                    "医療機関"                         "小計" 
                            si                          total 
"民間検査会社（主に自費検査）"                           "計" 

以降の処理のために date 列を関数 lubridate::date() で date 型に変換する． 列の変換・追加などには関数 dplyr::mutate() を用いる．

(pcr_data <- mutate(pcr_data, date = date(date)))

# A tibble: 1,175 × 10
   date        niid   ciq    hc    ai  univ    mi   sub    si total
   <date>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1 2020-02-18   472    75   398     0    79    NA  1024    NA  1024
 2 2020-02-19    15    68   609     0     0    NA   692    NA   692
 3 2020-02-20    20    15   758     0     0    NA   793    NA   793
 4 2020-02-21   261   188   902   132   108    NA  1591    NA  1591
 5 2020-02-22   341   127   677     2    19    NA  1166    NA  1166
 6 2020-02-23    53    72   529     0     0    NA   654    NA   654
 7 2020-02-24    22   103   460    17     0    NA   602    NA   602
 8 2020-02-25   195    38   695     0     0    NA   928    NA   928
 9 2020-02-26   267    19   970   149     0    NA  1405    NA  1405
10 2020-02-27   237    61  1101     0     0    NA  1399    NA  1399
# ℹ 1,165 more rows

Tip

以下，この項で用いた関数に関する補足．

#' 関数 print() を用いると表示する行数を指定できる
print(pcr_data, n = 5) # 全ては n = Inf
#' 日本語を含むファイルでは文字化けが起こった場合は以下で対応する
#' 関数 readr::guess_encoding() でファイルの文字コードを推測する
guess_encoding("data/pcr_case_daily.csv")
#' "UTF-8" であると 1 の信頼度で認識される
#' 文字コードを指定して読み込む場合は以下のように記述する
pcr_data <- # 文字コードとして UTF-8 を指定
    read_csv(file = "data/pcr_case_daily.csv",
             locale = locale(encoding = "utf-8"))
#' その他の文字コードとしては "sjis", "shift-jis", "shift_jis", 
#' "cp932"(拡張文字を含む)などを大文字小文字は区別せず指定できる
#' URLを指定して読み込むこともできる 
pcr_data <- # 更新される情報を追跡する場合に利用を推奨
    read_csv("https://www.mhlw.go.jp/content/pcr_case_daily.csv")
#' 列名の変更にはいろいろな方法があるので適宜使用する
#' 読み込み時に行う方法
pcr_data <- read_csv("data/pcr_case_daily.csv",
                     skip = 1, # 列名の行を読み飛ばす
                     col_names = c("date","niid","ciq","hc","ai",
                                   "univ","mi","sub","si","total"))
pcr_data <- mutate(pcr_data, date = date(date)) # date型に変更
#' 関数 dplyr::rename() を使う方法
pcr_data <- read_csv("data/pcr_case_daily.csv") # そのまま読み込む
(pcr_colnames <- set_names(names(pcr_data), # 新旧の列名に対応するベクトルを作成
                           c("date","niid","ciq","hc","ai",
                             "univ","mi","sub","si","total")))
pcr_data <- rename(pcr_data, all_of(pcr_colnames)) # 列名を変更
pcr_data <- mutate(pcr_data, date = date(date)) # date型に変更

データフレームの操作

問題

pcr_case_daily.csv から以下の条件を満たすデータを取り出しなさい．

1. 医療機関 (mi) での検査件数が2000を越えたときの国立感染症研究所 (niid) と医療機関 (mi) のデータ．
2. 大学等 (univ) と医療機関 (mi) でともに検査件数が2000を越えたデータ．
3. 2020年3月の各機関(sub,total は集計なので除く)の検査件数データ．

解答例

それぞれ以下のようにすれば条件に合致したデータを抽出することができる．

pcr_data |>            # データフレーム
    filter(mi > 2000) |> # 行の条件による絞り込み
    select(c(niid,mi))   # 列の選択 select(niid,mi) としても良い

# A tibble: 977 × 2
    niid    mi
   <dbl> <dbl>
 1     0  3510
 2     0  2636
 3     0  3459
 4     0  3694
 5     0  3437
 6    83  2579
 7     0  4247
 8    36  3469
 9     0  3725
10    20  3643
# ℹ 967 more rows

pcr_data |>
    filter(univ > 2000 & mi > 2000) # 複合的な条件の指定

# A tibble: 746 × 10
   date        niid   ciq    hc    ai  univ    mi   sub    si total
   <date>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1 2020-08-12     0   107  5805 15264  2152  3694 27022    NA 27022
 2 2020-08-17     0    11  3193 14606  2303  4247 24360    NA 24360
 3 2020-08-18    36     4  5631 15884  2322  3469 27346    NA 27346
 4 2020-08-19     0    77  5505 16463  2412  3725 28182    NA 28182
 5 2020-08-20    20    37  4940 17632  2397  3643 28669    NA 28669
 6 2020-08-24     0     3  3123  6599  2482  3454 15661    NA 15661
 7 2020-08-25     0     0  4740 14602  2291  3512 25145    NA 25145
 8 2020-08-26     0   111  4611 16481  2496  3488 27187    NA 27187
 9 2020-08-27    27     0  5183 16751  2071  3412 27444    NA 27444
10 2020-08-31     0     0  2948 16133  2714  4578 26373    NA 26373
# ℹ 736 more rows

pcr_data |>
    filter(date >= "2020-03-01" & date < "2020-04-01") |> # 日付の範囲での指定
    select(!c(sub,total)) # 列の削除は select(-c(sub,total)), select(-sub,-total) も可

# A tibble: 31 × 8
   date        niid   ciq    hc    ai  univ    mi    si
   <date>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1 2020-03-01     0    22   544   116     0    NA    NA
 2 2020-03-02   125    24  1038    15     0    NA    NA
 3 2020-03-03    20    27  1718    11     0    NA    NA
 4 2020-03-04    71    12  1398     7     0    NA    NA
 5 2020-03-05    83    11  1579     5     0    NA    NA
 6 2020-03-06   107    35  1786     5    62    68    NA
 7 2020-03-07    18    13  1616    11    20    44    NA
 8 2020-03-08     0    10   579     2     8     9    NA
 9 2020-03-09     0    27  1072     4    97    48    NA
10 2020-03-10     8    57  1763    97    65    62    NA
# ℹ 21 more rows

データフレームの集約

問題

pcr_case_daily.csv について以下の集計を行いなさい．

1. 各機関でのPCR検査件数の最大値．
2. 2021年の各機関でのPCR検査件数の月ごとの最大値．

ヘルプを用いて datasets::mtcars の内容を調べた上で以下の集計を行いなさい．

1. 気筒数 (cyl) ごとに排気量 (disp) の最大値，最小値．
2. 気筒数 (cyl) とギア数 (gear) ごとの燃費 (mpg) の平均値．

解答例

pcr_case_daily.csv において各機関でのPCR検査件数の最大値を求めるには基本的には以下のようにすればよい．

pcr_data |> summarise(across(!date, max))

# A tibble: 1 × 9
   niid   ciq    hc     ai  univ    mi    sub    si total
  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>
1   517  1733 16973 161077 12731    NA 272018    NA    NA

NAが含まれる列については最大値求めることができないので，最大値の計算で NA (欠損データ) を除く必要がある．

pcr_data |> # max の無名関数を利用する方法
    summarise(across(!date, \(x) max(x, na.rm = TRUE)))

# A tibble: 1 × 9
   niid   ciq    hc     ai  univ     mi    sub     si  total
  <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
1   517  1733 16973 161077 12731 136570 272018 131822 371225

pcr_data |> # package::purrr の lambda 式を利用する方法
    summarise(across(!date, ~ max(.x, na.rm = TRUE)))

# A tibble: 1 × 9
   niid   ciq    hc     ai  univ     mi    sub     si  total
  <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
1   517  1733 16973 161077 12731 136570 272018 131822 371225

2021年の月ごとの各機関でのPCR検査件数の最大値は以下のようになる．

pcr_data |>
    filter(year(date) == 2021) |>
    group_by(month(date)) |>
    summarise(across(!date, max))

# A tibble: 12 × 10
   `month(date)`  niid   ciq    hc     ai  univ    mi    sub     si  total
           <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>  <dbl>  <dbl>  <dbl>
 1             1     0    23 12042  58398  6556 23417  91261     NA  91261
 2             2     0    38  7032  39353  5533 22878  68847     NA  68847
 3             3     0     3  5117  66274  5005 23983  95939     NA  95939
 4             4     1     4  8866  75440  5121 30378 115724     NA 115724
 5             5     0    20 11181  89381  6392 35379 130410     NA 130410
 6             6     0    28  5433  68911  5439 28268  95330     NA  95330
 7             7    10    21  7739  63796  6500 36674 104740     NA 174233
 8             8     1    23 13200 109734  7397 54257 168857 110172 272564
 9             9     0    20 11602  94414  7472 43944 150671  90013 239984
10            10     0    21  2294  56724  5521 29978  86137  76467 159101
11            11     0    21  1022  31346  5706 27344  55361  71597 123910
12            12     0     6  2452  24376  5733 26647  51718  66888 112438

datasets::mtcars の集計は以下のとおりである．

mtcars |>
    group_by(cyl) |>
    summarise(max_disp = max(disp))

# A tibble: 3 × 2
    cyl max_disp
  <dbl>    <dbl>
1     4     147.
2     6     258 
3     8     472 

mtcars |>
    group_by(cyl) |>
    summarise(min_disp = min(disp))

# A tibble: 3 × 2
    cyl min_disp
  <dbl>    <dbl>
1     4     71.1
2     6    145  
3     8    276. 

mtcars |> # まとめて計算することも可能
    group_by(cyl) |> # 列名の作られ方に注意
    summarise(across(disp, list(max = max, min = min)))

# A tibble: 3 × 3
    cyl disp_max disp_min
  <dbl>    <dbl>    <dbl>
1     4     147.     71.1
2     6     258     145  
3     8     472     276. 

mtcars |>
    group_by(cyl, gear) |>
    summarise(mpg = mean(mpg))

# A tibble: 8 × 3
# Groups:   cyl [3]
    cyl  gear   mpg
  <dbl> <dbl> <dbl>
1     4     3  21.5
2     4     4  26.9
3     4     5  28.2
4     6     3  19.8
5     6     4  19.8
6     6     5  19.7
7     8     3  15.0
8     8     5  15.4

Tip

グループは既定値では順次解除されるので以下のような集計も可能である．

mtcars |> 
    group_by(cyl, gear) |>
    summarise(mpg = mean(mpg)) |> # cyl のグループは残っている
    summarise(mpg = max(mpg)) # cyl ごとに平均値の最大値を求める

基本的なグラフの描画

問題

pcr_case_daily.csv を用いて以下の描画を行いなさい．

1. 検疫所 (b)，地方衛生研究所.保健所 (c)，民間検査会社 (d) における検査件数の推移．
2. 民間検査会社 (d)，大学等 (e)，医療機関 (f) での検査件数の関係 (散布図)．

解答例

書き方はいろいろあるので，以下はあくまで一例である．

if(Sys.info()["sysname"] == "Darwin") { # MacOSの場合には日本語フォントを指定する
    theme_update(text = element_text(family = "HiraginoSans-W4"))}
pcr_data |>
    select(c(date,ciq,hc,ai)) |> # 描画対象の列を抽出
    pivot_longer(!date, names_to = "organ", values_to = "nums") |> # 
    ggplot(aes(x = date, y = nums, colour = organ)) +
    geom_line() +
    labs(x = "日付", y = "検査件数")

y軸を対数表示にする場合は以下のようにすればよい．

pcr_data |>
    select(c(date,ciq,hc,ai)) |> # 描画対象の列を抽出
    pivot_longer(!date, names_to = "organ", values_to = "nums") |> # 
    ggplot(aes(x = date, y = nums, colour = organ)) +
    geom_line() +
    scale_y_log10() + # y軸を対数表示 (log10(0)=-Inf の警告が出る場合がある)
    labs(x = "日付", y = "検査件数")

散布図を並べるために GGally パッケージを利用する．

pcr_data |>
    select(c(ai,univ,mi)) |> # 描画対象の列を抽出
    GGally::ggpairs(columnLabels = pcr_colnames[c("ai","univ","mi")]) # ラベルを渡す

擬似乱数

問題

ヘルプを用いて以下の関数を調べよ．

1. 関数 base::sample()
2. 関数 stats::rbinom()
3. 関数 stats::runif()
4. 関数 stats::rnorm()
5. 関数 base::set.seed()

以下の試行を実装してみよ．

1. サイコロを10回振る．
2. 4枚のコインを投げたときの表の枚数．

解答例

それぞれの関数の基本的な使い方は以下の例のようになる(問題の試行も含まれる)．

#' 関数sampleの使い方
(x <- 1:10)   # サンプリング対象の集合を定義

 [1]  1  2  3  4  5  6  7  8  9 10

set.seed(123) # 乱数のシード値(任意に決めてよい)を指定
sample(x, 5)                 # xから5つの要素を重複なしでランダムに抽出

[1]  3 10  2  8  6

sample(x, 5, replace = TRUE) # xから5つの要素を重複ありでランダムに抽出

[1]  5  4  6  9 10

sample(x, length(x))         # xの要素のランダムな並べ替え

 [1]  5  3  9  1  4  7 10  6  8  2

sample(1:6, 10, replace = TRUE)             # サイコロを10回振る実験の再現

 [1] 2 2 1 6 3 4 6 1 3 5

sample(1:6, 10, prob = 6:1, replace = TRUE) # 出る目の確率に偏りがある場合

 [1] 1 1 2 2 2 1 1 1 2 1

#' 関数rbinomの使い方
rbinom(10, size = 4, prob = 0.5) # 表(1)の出る確率が0.5にコインを4枚投げる試行を10回

 [1] 3 0 2 3 1 2 1 1 3 3

rbinom(20, size = 4, prob = 0.2) # 個数を20, 確率を0.2に変更

 [1] 0 1 0 0 0 1 1 1 1 1 1 1 1 1 0 1 0 0 1 0

#' 関数runifの使い方
runif(5, min = -1, max = 2) # 区間(-1,2)上の一様乱数を5個発生

[1] -0.6665937 -0.2691416  1.0041668  0.2529403  1.3645875

runif(5)                    # 指定しない場合は区間(0,1)が既定値

[1] 0.1028646 0.4348927 0.9849570 0.8930511 0.8864691

#' 関数rnormの使い方
rnorm(10, mean = 5, sd = 3) # 平均5，分散3^2の正規乱数を10個発生

 [1]  2.196845  6.181126  6.210894  2.340690  1.043187  5.086532  3.703611
 [8] 10.069618  8.685178  5.828070

rnorm(10)                   # 指定しない場合は mu=0, sd=1 が既定値

 [1] -1.0489755 -0.5208693  1.6232025 -1.0700682  1.6858872 -0.2416898
 [7] -0.4682005 -0.7729782  2.1499193 -1.3343536

#' 関数set.seedについて
set.seed(1) # 乱数の初期値をseed=1で指定
runif(5) 

[1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819

set.seed(2) # 乱数の初期値をseed=2で指定
runif(5)    # seed=1の場合と異なる結果

[1] 0.1848823 0.7023740 0.5733263 0.1680519 0.9438393

set.seed(1) # 乱数の初期値をseed=1で指定
runif(5)    # 初めのseed=1の場合と同じ結果

[1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819

双六ゲーム

問題

以下の簡単な双六ゲームを考える．

• ゴールまでのます目は100とする．
• さいころを振り出た目の数だけ進む．
• ゴールに辿り着くまで繰り返す．

さいころを振る回数の分布がどうなるか実験を行いなさい．

解答例

双六の試行を行う関数を作成する．

mc_trial <- function(){
    step <- 0 # 最初の位置
    num <- 0  # さいころを振る回数
    while(TRUE){ # 永久に回るループ
        step <- step + sample(1:6, 1) # さいころを振る
        num <- num + 1 # 回数を記録
        if(step >= 100) { # ゴールしたか?
            return(num) # 回数を出力して関数を終了
        }
    }
}

Tip

同じ試行でも関数の作り方はいろいろある． 例えば100回サイコロを振れば必ずどこかで100を越えるので，計算は無駄があるが条件分岐のない関数を考えることもできる．

mc_trial <- function() {
    which.max(cumsum(sample(1:6, 100, replace = TRUE)) >= 100)
} # 関数 which.max() は初めて TRUE(1) になった場所を返す

試しに10回の試行を実行してみる．

for(i in 1:10) print(mc_trial())

[1] 29
[1] 27
[1] 27
[1] 31
[1] 29
[1] 33
[1] 30
[1] 26
[1] 27
[1] 31

確率シミュレーション(Monte-Carlo実験)を行う．

set.seed(12345) # 必要に応じて乱数のシードを設定する
mc_num <- 10000 # 実験回数を設定 
mc_data <- replicate(mc_num, mc_trial()) 
summary(mc_data) # 簡単な集計

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  20.00   27.00   29.00   29.09   31.00   41.00 

ヒストグラムを作成して結果を視覚化する．

tibble(x = mc_data) |> # ヒストグラムを出力
    ggplot(aes(x)) + 
    geom_histogram(binwidth = 1,
                   fill = "slateblue", alpha = 0.5, # 塗り潰しの色
                   colour = "slateblue")  # 縁の色