準備

まず、今回利用するパッケージを読み込む。インストール済みでない場合は、読み込む前にインストールする。

# 必要ならまずインストール
# install.packages('dplyr', dependencies = TRUE) 
# install.packages('ggplot2', dependencies = TRUE)
library("dplyr")
## 
## Attaching package: 'dplyr'
## 
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## 
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library("ggplot2")
# 次の行はMac用（Windows上では実行しないように）
theme_set(theme_gray(base_size = 12, base_family = "HiraKakuProN-W3"))

地球儀問題の検討

講義で紹介した地球儀問題 (McElreath 2016: Ch.2) を検討しよう。

地球儀のトスを1回行う度にベイズ更新を行う場合

事前の妥当性

まず、事前の妥当性を図示する。すべての \(\theta\)（水の割合）が等しく妥当だと仮定する。

# 母数 theta を含むデータフレームを用意する
delta <- .001
globe <- data.frame(theta = seq(0, 1, by = delta))
# すべての theta の事前の妥当性 (prior) を1に設定する
globe <- mutate(globe, prior = rep(1, n()))
# このデータセットの中身を確認したい場合は、次の2行を実行（コメント'#'を外す）
# head(globe)
# tail(globe)
# 図を表示する
plt_prior <- ggplot(globe, aes(x = theta, y = prior)) +
    geom_line(color = "tomato") +
    xlab('Wの割合')　+ ylab('妥当性')
    ylim(0, 2.5)  # 縦軸の範囲を指定
## <ScaleContinuousPosition>
##  Range:  
##  Limits:    0 --  2.5
print(plt_prior)

ベイズ更新

まず、1回目の観測値 W でベイズ更新を行う。 1回のトスでWが出る確率は Wの割合 \(\theta\) なので、尤度は\(\theta\) である。 Rで尤度を設定し、事後確率を計算する。

globe$lik <- globe$theta
globe$post <- with(globe, lik * prior)
globe$post <- globe$post / sum(delta * globe$post)

この結果を図示する。

plt_post <- ggplot(globe, aes(x = theta)) +
    geom_line(aes(y = prior)) +
    geom_line(aes(y = post), color = "red", size = 1.5) +
    xlab('Wの割合') + ylab('妥当性')
plt_post

この図では、黒い線がデータを観察する前の妥当性（事前確率）を、赤い線がデータ観察後の妥当性（事後確率）を表している。

2回目以降の更新も同様に行う。各トスにおける尤度（データ生成過程は）、次の数式で表現できる。 \[L(\theta | w) = \Pr(w|\theta) = \theta^w (1 - \theta)^{1-w}.\] ただし、観察された結果がWのとき\(w=1\)、Lのとき\(w=0\)とする。 つまり、Wがでたときの尤度は\(\theta\)、Lが出たときの尤度は \(1-\theta\) である。 これは、ベルヌーイ分布 (Bernoulli distribution) の確率（確率質量）である。

ベイズ更新は新しい観測値を得るたびに同じ作業を繰り返し行うので、作業内容をまとめた関数を作ろう（2回以上同じことを繰り返すときは、関数を使ったほうが楽である）。

R で自前の関数を作るときは、function() を使う。 ここでは、ベイズ更新の過程を globe_update という名前の関数として定義してみよう。

関数を作るときにまず必要なのは、引数（ひきすう、argument）である。 これは、関数の入力値となるもので、ここでは observation、 data_frame、first の の3つの引数を使うことにする。この関数を利用する際には、これらの引数を入力する必要がある。ただし、3つ目の first のみ first = FALSE としてデフォルトの値を指定しておく。関数を利用する際に first を指定しないと、デフォルトの値である FALSE が利用される。

次に、関数が最終的に何を出力するか決める必要がある。 関数の出力のことを、返り値または戻り値 (return) と呼び、これをreturn() で指定する。ここでは、入力されたデータフレームに更新された事後確率を新たな変数として加えたデータフレームを返すことにする。関数は、return() に到達すると終了する。したがって、関数の定義が10行書かれていても、第1行にあるreturn()が評価されれば、残りの行は無視される。

指定された引数を使い、私たちが必要とする返り値を得るための関数（の1例） は次の通りである。現時点でコードを完全に理解する必要はないが、1行ずつ実行してみれば、各行の内容はそれほど難しくないだろう。

globe_update <- function(observation, data_frame, first = FALSE) {
    # 地球儀問題でベイズ更新を行うための関数
    # 引数: observation = {W, L}
    #       data_frame = thetaベクトルと事前確率を含むデータフレーム
    #       first = {TRUE (1回目の更新), FALSE (2回目以降の更新)}
    # 返り値: df = 事後確率を加えたデータフレーム

    df <- data_frame
    if (!first) { # 2回目以降の更新：前回の事後確率を事前確率として使う
        df$prior <- df$post
    }
    prior <- df$prior
    if (observation == "W") {
        lik <- df$theta
    }
    else {
        lik <- 1 - df$theta
    }
    post <- lik * prior
    delta <- 1 / (nrow(df) - 1)
    df$post <- post / sum(delta * post)
    return(df)
}

これで、関数が定義できた。

関数について現時点で注意すべき点として、関数内での定義された変数には関数の外からアクセスできないという点がある。それらの変数は、あくまで関数内でのみ利用されるものである（この点について気になるなら、「変数のスコープ」をキーワードにしてネットで検索してみよう）。

この関数を使い、再び1度目の更新から実行してみよう。

# 母数 theta を持つデータフレームを作る
globe_0 <- data.frame(theta = seq(0, 1, by = .001))
# すべての theta の事前確率 (prior) を1に設定する
globe_0$prior <- rep(1, nrow(globe_0))
# 地球儀問題のために作った関数を使ってベイズ更新！
globe_1 <- globe_update(observation = "W", data_frame = globe_0, first = TRUE)

図は、上で作った図のデータフレームのみアップデートする。

plt_post %+% globe_1

このように、ggplot で1度作ったものと同様の図を、異なるデータフレームを使って作り直すときは、もう1度すべてのコードを書く必要はなく、%+% を使ってアップデートすればよい。

続いて、2回目の更新を行う。観察されたのはLである。1回目の更新ではないので、first 引数はデフォルト値を使う。関数の定義と同じ順番で引数を並べるときは、引数名を省力できる。

globe_2 <- globe_update("L", globe_1)
plt_post %+% globe_2

3回目以降も同様に更新する。 3回目はW。引数名を明示すれば、順番は変えてもよい。

globe_3 <- globe_update(data_frame = globe_2, observation = "W")
plt_post %+% globe_3

4回目はW。

globe_4 <- globe_update("W", globe_3)
plt_post %+% globe_4

5回目はW。

globe_5 <- globe_update("W", globe_4)
plt_post %+% globe_5

6回目はL。

globe_6 <- globe_update("L", globe_5)
plt_post %+% globe_6

7回目はW。

globe_7 <- globe_update("W", globe_6)
plt_post %+% globe_7

8回目はL。

globe_8 <- globe_update("L", globe_7)
plt_post %+% globe_8

9回目はW。

globe_9 <- globe_update("W", globe_8)
plt_post %+% globe_9

更新の順番は影響するか？

上のデータは、\(D_1 = (W, L, W, W, W, L, W, L, W)\)であった。 もし、データが\(D_2 = (L, L, L, W, W, W, W, W, W)\)だったら、どのような更新になるだろうか。また、最終的に得られる妥当性は異なるものになるだろうか。各自で確かめてみよう！

また、観測した値を一括してベイズ更新を行うと、どのような結果になるだろうか。これについては以下のグリッド近似を参照されたい。

グリッド近似による事後確率の計算

グリッド近似では、母数（パラメタ）の取り得る範囲内からいくつかの点（グリッド）を選び、選ばれた各点で事後確率を計算する（上で実行したのも実はグリッド近似である）。具体的には、以下のステップが必要である。

  1. グリッドを定義する。母数の範囲内でいくつの点（グリッド）をとるか決め、それぞれの点で母数がとる値を定義するベクトル (vector) を決める。
  2. 各点における母数の事前確率を決める（計算する）。
  3. 各点で尤度を計算する。
  4. 各点で尤度と事前確率を掛け合わせる（標準化前の事後分布を求める）。
  5. 4 で得られた値を標準化し、各点の事後確率を計算する。

地球儀問題を例にして、Rでグリッド近似を行ってみよう。

1. グリッドの定義

まず、グリッドを定義する。推定する母数\(\theta\) は0から1の間の実数をとる。この範囲に、5つの点を等間隔に取ってみよう。 そのために、seq() 関数を使う。

n_grid <- 5  # 点の数を指定する
theta_grid <- seq(from = 0, to = 1, length.out = n_grid)

これで、0から1までの間に5点とることができた。 念のため、グリッド中身を確認してみよう。

theta_grid
## [1] 0.00 0.25 0.50 0.75 1.00

このように、5点が等間隔で取られている。

2. 事前確率の設定

ここでは、どの値の妥当性も、データを観察する前は同じであると仮定しよう。 そうすると、グリッド上の各母数の値を定数で表せばよいので、rep() 関数を使って、グリッドの数と同じ数の1が並ぶベクトルを作る。

( prior <- rep(1, n_grid) )
## [1] 1 1 1 1 1

（このように、変数を定義する （<- で値を割りあてる）とき、全体をカッコで括ると、割り当てと同時中身が表示される 。）

3. 尤度の計算

地球儀問題のデータは、二項分布 (binomial distribution) によって生み出されていると考えることができる。 \(D=(W, L, W, W, W, L, W, L, W)\) というデータは、9回のトスのうちWが6回、Lが3回出たことを示す。また、母数 \(\theta\) は W の割合を表しているので、尤度は、 \[L(\theta | w=6, n=9) = \Pr(w=6 | n=9, \theta) = \binom{9}{6}\theta^6 (1 - \theta)^3\] である。

練習のために、この数式を表す自前の尤度関数を作ってみよう。

my_lik <- function(w, n, theta) {
    # 引数: w = W が出た回数
    #       n = トスを行った回数
    #       theta = 各トスでWが出る確率
    # 返り値: lik = 尤度
    lik <- choose(n, w) * theta^w * (1 - theta)^(n - w)
    return(lik)
}

ここでは、w, n, theta の3つを引数（ひきすう, argument）に指定し、関数にmy_lik という名前をつけた。 この関数は与えられた引数を使って定められた計算をし、lik を私たちに返す。

この関数を使って尤度を計算してみよう。

lik_0 <- my_lik(w = 6, n = 9, theta = theta_grid)
lik_0
## [1] 0.000000000 0.008651733 0.164062500 0.233596802 0.000000000

このように、各点における尤度を計算することができる。

ここで定義した尤度は二項分布であり、Rには二項分布のためのdbinom() という関数が用意されているので、その関数を利用してみよう。 この関数の引数は、x = 「成功」の回数（地球儀問題ではWの数）、n = 試行の回数（同じくトスの回数）、prob=1回の試行における成功確率（地球儀表面の水の割合）の3つである。

lik <- dbinom(6, size = 9, prob = theta_grid)
lik
## [1] 0.000000000 0.008651733 0.164062500 0.233596802 0.000000000

このように、より簡単に尤度が求められる。自分で定義した尤度関数で計算したものと各値が同じである。

rbind(lik_0, lik)
##       [,1]        [,2]      [,3]      [,4] [,5]
## lik_0    0 0.008651733 0.1640625 0.2335968    0
## lik      0 0.008651733 0.1640625 0.2335968    0

4と5. 事後確率の計算

事後確率を求めるために、まず尤度と事後分布をかける。

post <- lik * prior
post
## [1] 0.000000000 0.008651733 0.164062500 0.233596802 0.000000000
sum(post)
## [1] 0.406311

この時点では、事後確率が標準化されていない。つまり、事後確率の合計が1にならない。

事後「確率」にするためには、合計が1になるように調整すれば良いので、sum() で求められる合計で各値を割る。

post <- post / sum(post)
post
## [1] 0.00000000 0.02129338 0.40378549 0.57492114 0.00000000
sum(post)
## [1] 1

これで事後確率を求めることができた。

最後に、事後確率を図示しよう。 ggplot で作図するため、必要な変数を含んだデータフレームをdata.frame() で作る。 私たちが作りたい図は、縦軸に事後確率、横軸にWの割合を示すものなので、theta_grid とpost を変数に持つデータフレームを作り、df_1 という名前をつける。変数名は、それぞれtheta とposterior に変更する。

df_1 <- data.frame(theta = theta_grid, posterior = post)
df_1
##   theta  posterior
## 1  0.00 0.00000000
## 2  0.25 0.02129338
## 3  0.50 0.40378549
## 4  0.75 0.57492114
## 5  1.00 0.00000000

このデータフレームを使ってggplotで事後分布を図示しよう。

plt_grid_1 <- ggplot(df_1, aes(x = theta, y = posterior)) +
    geom_line(color = "gray") +  # 事後確率を線で繋ぐ
    geom_point(size = 4)         # グリッド上の事後確率を点で表示する
plt_grid_1

事後確率を自分で確認するにはこれで十分だが、ラベルも書き直してみよう。

plt_grid_1 <- plt_grid_1 +
    labs(x = "水面の割合", y = "事後確率", title = paste0(n_grid, "点グリッド"))
plt_grid_1

ここで使ったpaste0() 関数は、中身の文字列を隙間なく繋ぐために使う。間に何か文字（あるいはスペース）を入れたいときは、paste() を使う。例えば、

six <- 6
paste0(six, "-inch")
## [1] "6-inch"
paste(six, "inch", sep = "-")
## [1] "6-inch"
paste("There are", six, "students.", sep = " ")
## [1] "There are 6 students."

のように使う。

表示された図はレポートなど再利用するために、PDFファイルに保存する。 （他の形式でも構わないが、LaTeX での使用を想定すると、PDFが最も望ましい。） Mac では、以下のようにする。

quartz(file = "grid-approx-5points.pdf", type = "pdf", family = "sans", 
       height = 6, width = 6)
print(plt_grid_1)
dev.off()
## quartz_off_screen 
##                 2

Windows では次のようにする。

pdf(file = "grid-approx-5points.pdf", width = 6, height = 6)
print(plt_grid_1)
dev.off()

quartz() またはpdf() でグラフィックデバイスを開き、そこに図を書き込み、dev.off() でデバイスを閉じる。 dev.off() を忘れないように注意すること。図の幅 (width) と高さ (height) はインチで指定する。 file で指定した名前のファイルが、現在の作業ディレクトリに保存される。

R Studio では、図のすぐ上にある “Export” ボタンから保存することもできるが、コマンドで保存する方がおすすめ（その方が結果的に楽なので）である。

グリッドの数を変えてみる

同様のグリッド近似を、グリッドの数を変えて実行してみよう。 変える必要があるのは、n_grid だけである。

n_grid <- 20
theta_grid <- seq(from = 0, to = 1, length.out = n_grid)
prior <- rep(1, n_grid)
lik <- dbinom(6, size = 9, prob = theta_grid)
post <- lik * prior
post <- post / sum(post)

ここで、先ほどと同じように事後確率を図示してみよう。 上のコードを真似て図を作ると。

df_2 <- data.frame(theta = theta_grid, posterior = post)
plt_grid_2 <- ggplot(df_2, aes(x = theta, y = posterior)) +
    geom_line(color = "gray") +
    geom_point(size = 4) +
    labs(x = "水面の割合", y = "事後確率", title = paste0(n_grid, "点グリッド"))
plt_grid_2

これで図は作れるが、やや面倒である。先ほどの図 plt_grid_1 とこの plt_grid_2 を作ったコードを比較してみると、違うのは利用したデータフレームだけで、残りは全く同じである。 そこで %+% を用いて図に利用するデータフレームだけをアップデートする。

plt_grid_2 <- plt_grid_1 %+% df_2
plt_grid_2

この図は、タイトルが望んでいるものとは違う。こうなってしまった理由は、タイトルに使う n_grid がデータフレームの中に含まれておらず、n_grid の変更がアップデートに反映されていないためである。 そこで、次のようにタイトルを書き換える。

plt_grid_2 <- plt_grid_2 + labs(title = paste0(n_grid, "点グリッド"))
plt_grid_2

これで、20個の点を使ったグリッド近似ができた。

グリッドの数を変えて、グリッド近似の練習をしてみよう！



授業の内容に戻る