準備

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

# 必要ならまずインストール
# 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"))

dplyr パッケージについては、第9回「データの編集」で詳しく説明する。 基本的な使い方については、TA の宋さんが書いた解説 や cheat sheet (PDF) を参照されたい。

Rを使った [擬似] 乱数 ([pseudo-]random numbers) の生成

Rを使うと、様々な方法で乱数 (random numbers) 生成することができる。 ただし、完全にランダムに数ができるわけではないので、より正確には擬似乱数 (pseudorandom numbers) が生成できる。 今回はまず、基本的な乱数の作成方法を学習しよう。

特定のベクトルからの無作為抽出

あるベクトルに含まれる値を無作為に選びたいときは、sample() を使う。 この関数で必ず指定する必要がある引数は、抽出の元となるベクトル x と、そこからいくつの値を抽出するかを指定する size の2つである。

まず、1から10までの整数ベクトルx を用意し、そこから2つの値を無作為に抽出してみよう。

(x <- 1:10) 
##  [1]  1  2  3  4  5  6  7  8  9 10
sample(x, size = 2)
## [1] 6 5

このように、1から10の整数の中から2つの値が取り出される。

取り出された値が本当に無作為に選ばれているか確かめるために、 何度か同じコードを実行してみよう。

sample(x, size = 2)
## [1] 2 9
sample(x, size = 2)
## [1] 10  5
sample(x, size = 2)
## [1] 9 4
sample(x, size = 2)
## [1] 9 3
sample(x, size = 2)
## [1] 10  2

このように、毎回異なるペアが抽出される。もちろん、たまたままったく同じ値のペアが抽出されることもある。

無作為に選ぶことが本来の目的なので、毎回異なる値が出るのは望ましいことである。 しかし、レポートや論文などを書く際、数字が毎回変わってしまうと都合が悪いことがある。 例えば、乱数を図にして説明する際、図を作り直すたびに数字が変わってしまうと、その度に本文の数値例も書き換える必要がある。 （ただし、R Markdown を使えば、このような心配をする必要はあまりない。詳しくは、R Markdown の使い方の回で。） そのような自体を避けるため、乱数を生み出す種 (seed) を指定することで毎回同じ値を得ることができる。 種の指定には、set.seed() を使う。 種の指定の仕方に特に決まりはないので、適当に数値を与える。ただし、毎回同じ数字を使うと、乱数を使っているとは言えなくなるので注意されたい。

set.seed(12345)  # 乱数の種を指定
sample(x, size = 2)
## [1]  8 10
set.seed(12345)  # 上と同じ種を指定
sample(x, size = 2)
## [1]  8 10
set.seed(12345)  # 上と同じ種を指定
sample(x, size = 2)
## [1]  8 10

このように、何度でも同じ結果が出せる。

抽出の元となるベクトルの値は、数字でなくてもよい。例えば、赤い玉 (R) と白い玉 (W) から無作為に1つを選びたいなら、 次のようにすればよい。

balls <- c('R', 'W')
sample(balls, size = 1)
## [1] "W"

ボールは2つしかないので、2つより多い数を選ぼうとするとエラーが出る。

sample(balls, size = 10)
## Error in sample.int(length(x), size, replace, prob): cannot take a sample larger than the population when 'replace = FALSE'

2つのボールが入った袋から1つを無作為に取り出し、そのボールを袋に戻してからもう1度取り出すという作業を繰り返すなら、2つより多くボールを取り出すことができる。このように、1度取り出した値を元に戻してから次の抽出を行う方法を、復元抽出 (sampling with replacement) と呼ぶ。sample() 関数では、replace = TRUE と指定することにより、復元抽出を実行することができる。 （すなわち、sample() はデフォルトで replace = FALSE となっており、非復元抽出 (sampling without replacement) を実行する。）

sample(balls, size = 10, replace = TRUE)
##  [1] "W" "R" "R" "R" "W" "W" "W" "R" "R" "W"

このように、2つしかボールがなくても、復元抽出により10個の観測値を得ることができる。

赤玉が4つ、白玉が1つ入っているバッグから復元抽出で無作為にボールを10回取り出す作業は、以下のように実行できる。

bag <- c('R', 'R', 'R', 'R', 'W')
sample(bag, size = 10, replace = TRUE)
##  [1] "R" "R" "R" "R" "R" "R" "W" "R" "R" "W"

この例では、赤と色の個数（あるいはその比）が正しく反映されたベクトル bag を用意したが、 各値がでる確率（比率）は prob で指定することもできる。prob の値は、デフォルトではすべての値が等確率で出るように設定されている。

赤と白が1つずつしか含まれていないベクトル balls から、赤と白の比率が4:1 で出るような無作為復元抽出を行うには、以下のようにする。

set.seed(2016-05-05)  # 下のコードと比較するために種を指定
sample(balls, size = 10, replace = TRUE, prob = c(4/5, 1/5))
##  [1] "W" "W" "R" "R" "R" "R" "R" "R" "R" "R"

このように、prob にベクトルの各値が1回の抽出で選ばれる確率を指定する。 このコードは、balls の第1要素である ‘R’ が出る確率が4/5、第2要素である ‘W’ が出る確率が1/5になるように書かれている。

prob に指定するのは確率でなくてもよい。各値の出やすさを重みで表すこともでできる。 赤と白の比を4:1 にするために、以下のようにしてもよい。

set.seed(2016-05-05)  # 上のコードと比較するために同じ種を指定
sample(balls, size = 10, replace = TRUE, prob = c(4, 1))    
##  [1] "W" "W" "R" "R" "R" "R" "R" "R" "R" "R"

先ほどと同じ結果が得られる。

ここで、私たちが袋の中に赤と白のボールがいくつずつ入っているか知らないとしよう。 そんなときは、袋からボールを無作為に復元抽出する作業を何度も繰り返せば、袋に含まれる赤いボールの比率を推定することができる。 具体的には、

  1. ボールを n 回取り出す
  2. 赤が出た回数を n で割る

という作業を行えばよい。実際に袋とボールを用意してこの作業を行うことも可能である。nが100程度であればやってもよいかもしれない。しかし、n > 1000 となるとかなり面倒だろう。そのようなときは、R で作業を行えばよい。コンピュータは退屈な繰り返し作業を任せるのに最適である。

試しに、n = 10,000 として赤いボールの比率を推定してみよう。

n <- 10000
experiment <- sample(bag, size = n, replace = TRUE) 
sum(experiment == 'R') / n
## [1] 0.7995

上のコードは、まず10,000回の抽出結果を experiment という名前のベクトルに保存する。したがって、experiment の中身は、‘R’ または ‘W’ で、要素の数は全部で 10,000 である。

その次のコードのうち、experiment == 'R' という部分で、experiment の各値が R に等しいかどうかを確かめている。このように、値が等しいかどうかを調べるときは、二重等号 == を使う。 結果として、等しければ TRUE が、等しくなければ FALSE が返される。つまり、experiment == 'R' の結果は、要素の数が10,000 で、各要素は TRUE または FALSE のベクトルである。 次に、そのベクトルの合計を sum() で計算している。計算に使われるとき、TRUE = 1、FALSE = 0 とされる。したがって、 このsum() はTRUE の数、すなわち 10,000 回のうち何回赤いボールが観察されたかを数えている。 最後に、この数を総抽出数 n = 10,000 で割って、赤いボールの比率を推定している。 この推定値は、真実の値（今回私たちは偶然にも袋の中身を知っていた）である 0.8 に近い値であり、サンプル抽出による推定がうまくいったようである。

確率分布からの乱数生成

正規分布 (normal distributions)

正規分布からの無作為抽出は、rnorm() で実行する。 この関数の使い方は、Rで関数を利用する で説明してある。

平均 10, 標準偏差 2 の正規分布から10個の値を無作為に抽出し、抽出した値の平均と標準偏差を計算してみよう。

norm_1 <- rnorm(10, mean = 10, sd = 2)
norm_1
##  [1]  9.855811 11.815302  7.347727 11.130927  8.656476 11.289818  9.974155
##  [8]  8.002123  8.572534 10.679008
mean(norm_1)
## [1] 9.732388
sd(norm_1)
## [1] 1.522902

このように、平均値は抽出元の分布の平均に近いが、標準偏差にはだいぶずれがある。

抽出した値をヒストグラムにしてみよう。

hist(norm_1, freq = FALSE, col = 'gray')

このヒストグラムは正規分布に見えないだろう。

次に、抽出する値の個数を増やし、もう1度試してみよう。

norm_2 <- rnorm(10000, mean = 10, sd = 2)
mean_norm_2 <- round(mean(norm_2), 2)
sd_norm_2 <- round(sd(norm_2), 2)
hist(norm_2, freq = FALSE, col = 'skyblue')
text(x = c(15, 15), y = c(.175, .15),
     labels  = c(paste0('mean = ', mean_norm_2),
                 paste0('sd = ', sd_norm_2)))

このように、抽出する値の数を増えせば、正規分布に見えるようになる。 また、平均、標準偏差とも、抽出元の分布に近づく。

標準正規分布 (standard normal distribution) からの抽出は次のようにできる。

norm_4 <- rnorm(10000) # rnorm(10000, mean = 0, sd = 1) と同じ
mean_norm_4 <- round(mean(norm_4), 2)
sd_norm_4 <- round(sd(norm_4), 2)
hist(norm_4, freq = FALSE, col = 'skyblue', ylim = c(0, .4))
text(x = c(2.5, 2.5), y = c(.35, .3),
     labels = c(paste0('mean = ', mean_norm_4), paste0('sd = ', sd_norm_4)))


一様分布 (uniform distributions)

最小値\(a\)、最大値\(b\)の連続一様分布 (uniform distribution) からの無作為抽出は、runif(n, min = a, max = b) で行う。\(n\) は抽出する値の数である。

一様分布の確率密度は、 \[p(x | a, b) = \frac{1}{b - a}\] であり、期待値は、 \[\mathrm{E}(x) = \int_a^b p(x | a, b) x dx = \frac{a + b}{2}\] である。

最小値\(-5\)、最大値\(15\) の連続一様分布を例に、無作為抽出を行ってみよう。 この分布の確率密度は \(1/[15 - (-5)] = 1/20 = 0.05\)、平均は \((-5 + 15)/2=5\) である。

unif_1 <- runif(10000, min = -5, max = 15)
hist(unif_1, freq = FALSE, col = 'tomato')

mean(unif_1)
## [1] 5.065406

10,000個の値を無作為に抽出したところ、確率密度が0.05 の一様分布に近いサンプルが得られ、その値の平均は約5になった。

最小値と最大値を指定しないと、runif() はデフォルトで最小値0、最大値1の一様分布を利用する。 この分布の確率密度は1、期待値は0.5 である。

unif_2 <- runif(10000)
hist(unif_2, freq = FALSE, col = 'tomato')

mean(unif_2)
## [1] 0.49654


二項分布 (binomial distributions)

サイズ（ベルヌーイ試行の数）Nでそれぞれ成功確率が \(\theta\) (theta) の二項分布から\(n\)個の値を無作為に抽出するには、 rbinom(n, size = N, prob = theta) を使う。

二項分布の確率質量は \[\Pr(x | N, \theta) = \binom{N}{x} \theta^x (1 - \theta)^{N - x},\] 期待値は \[\mathrm{E}(x) = \sum_{k=0}^N \Pr(k | N, \theta) k = N\theta,\] 分散は \[\mathrm{Var}(x) = \mathrm{E}\left([x - \mathrm{E}(x)]^2\right) = N \theta (1 - \theta)\] である。

例として、成功確率 0.3、サイズ\(N=10\) の二項分布を考える。この二項分布の取り得る値は、\(0, 1, 2, \dots, 10\) の11種類（10種類ではない！）ある。11個の各値について、上の式で確率質量が求められる。 また、期待値は \(10 \cdots 0.3 = 3\)、分散は \(10 \cdot 0.3 (1 - 0.3) = 2.1\) となる。 この分布から10,000個の値を無作為に抽出する。

binom_1 <- rbinom(10000, size = 10, prob = .3)
mean(binom_1)
## [1] 3.016
var(binom_1)
## [1] 2.125157

このように、平均、分散は元の分布に近いものが得られた。

このサンプルを図示してみよう。

hist(binom_1, col = 'orange')

この図は少し変である。ヒストグラムなのに棒と棒の間が空いているが、0と1の間だけ詰まっている。 よく考えてみると、二項分布は離散分布 (discrete distribution) なので、ヒストグラムは二項分布を表すのに適していない。離散分布には棒グラフ (bar plot) を使った方がよい。

そこで、ggplot を使って棒グラフを図示する。

df <- data.frame(binom_1)
plt_binom <- ggplot(df, aes(x = binom_1)) + geom_bar()
print(plt_binom)

これで棒グラフができたが、さらに2点で改良を加える。 まず、横道の値をみると、2.5 や7.5 など、二項分布では取り得ない値が表示されてしまっている。 これは、ggplot が binom_1 の値が連続値であると誤解しているせいである。 そこで、as.factor(binom_1) とすることにより、binom_1 をカテゴリ変数として認識させる。

次に、縦軸を回数 (count) ではなく、確率にする。今回のように大きなサンプルにおいて、確率はその値の回数を全体の回数で割ることによって推定される。ggplot は各値の回数を ..count.. に保存しているので、これを使って確率を求め、それを縦軸に指定する。 （..count.. や..density.. など、ggplot が内部で計算して図に使用する変数は、2つのピリオドが前後につく。）

plt_binom_better <- ggplot(df, aes(x = as.factor(binom_1))) +
    geom_bar(aes(y = (..count..) / sum(..count..)), fill = 'darkgreen') +
    labs(x = 'binom_1', y = 'probability')
print(plt_binom_better)

これで、それぞれの値が出る確率（より正確には、このサンプルで各値が出た割合）を図示できた。

二項分布のサイズをN = 1 にすれば、ベルヌーイ分布 (Bernoulli distribution) からの無作為抽出ができる。

bern_1 <- rbinom(10000, size = 1, prob = .3)
df$bern <- bern_1
plt_bern <- ggplot(df, aes(x = as.factor(bern_1), y = (..count..) / sum(..count..))) +
    labs(x = 'bern_1', y = 'probability')
print(plt_bern + geom_bar())

棒が少し太すぎるので、width で幅を調整する。

print(plt_bern + geom_bar(width = .5, fill = 'navy'))

成功確率0.3 のベルヌーイ試行を10,000回繰り返した結果、成功の比率が約0.3、失敗の比率が約0.7 のサンプルが得られた。



事後分布からのサンプリング

再び、地球儀問題 (McElreath 2016) を考える。 地球儀表面の水の割合を \(\theta\)とし、\(\theta\) の事前確率は一様分布 U(0, 1) に従うものとする。 また、観察されたデータは「9回のトスのうち6回が水」である。

事後分布からサンプルを用いた推論

この問題の事後分布は、グリッド近似によって以下のように求められる。 グリッド近似を後で簡単に実行できるよう、関数を作る。

globe_grid <- function(w, n, n_grid = 1000){
    # 地球儀問題のグリッド近似を行う
    # 引数: w = W の回数
    #       n = トスの回数
    #       n_grid = グリッドの数, デフォルトは1000
    # 返り値: 各グリッドの母数の値と事後確率を変数に持つデータフレーム 
    theta_grid <- seq(0, 1, length.out = n_grid)
    prior <- rep(1, n_grid)
    lik <- dbinom(w, size = n, prob = theta_grid)
    post <- lik * prior
    post <- post / sum(post)
    df <- data.frame(theta = theta_grid,
                     posterior = post)
    return(df)
}

作った関数にデータを与えてグリッド近似を行う。

df <- globe_grid(w = 6, n = 9)

この結果を図示する。

grid_plt <- ggplot(df, aes(x = theta, y = posterior)) + geom_line(color = 'firebrick')
print(grid_plt)

この図は、（標準化されていない [面積の合計が1にならない]） 事後分布を示している。 この図から、水の割合についてどのような推論をすることができるだろうか？ 確率密度曲線を数学的に解析することも可能だが、 この事後分布から無作為に値をたくさん抽出し、そのサンプル内の頻度を使って推論を行うことにする。

まず、事後分布から100000個の値を無作為に抽出する。

n_sim <- 100000
post_sim <- sample(df$theta, size = n_sim, prob = df$posterior, replace = TRUE)

抽出した値をプロットしてみよう。縦軸は、事後分布から抽出された \(\theta\) （水の割合）の値、横軸は抽出された順番（したがって、横軸は重要ではない）とする。

df_sim <- data.frame(post_sim) %>%  # データフレームを作る
    mutate(id = 1:n())              # 抽出順を示す id を作る
plt <- ggplot(df_sim, aes(x = id, y = post_sim)) +
    geom_point(color = 'firebrick', alpha = 1/20) +
    ylim(0, 1) +
    labs(x = 'サンプル番号', y = expression(theta))
print(plt)

ここでは、アルファブレンディングという手法を使って図を作っている。 アルファブレンディングとは、半透明な図を作る方法で、alpha の値を指定することで1点あたりの色の濃さを変える。ここでは、alpha = 1/20 としているので、点が20個重なる部分が通常の色の濃さ（透過率 = 0）になり、それ以下の点しかないところが半透明になる。 この方法は、プロットする点が多いときに有効である。

出来上がった図を見ると、0.4 から0.8 付近のあたりの色が濃くなっており、\(\theta\) の値がこの付近に集中していることがわかる。したがって、\(\theta\) の値は0.4から0.8あたりの値を取りやすいと推論することができる。

次に、抽出した値をヒストグラムにしてみよう。

hist_globe <- ggplot(df_sim, aes(x = post_sim, y = ..density..)) +
    geom_histogram(binwidth = .05, color = 'black') +
    labs(x = expression(theta), y = '確率密度')
print(hist_globe)

抽出された \(\theta\) の値は上のヒストグラムが表すとおり分布している。

ここで、事後分布が特定の範囲の値をとる確率を考えてみよう。 例えば、\(\theta < 0.5\) となる確率はいくつだろうか。 求めたい確率を図示すると次のようになる。

area_1 <- ggplot(df, aes(x = theta, y = posterior)) +
    geom_line(color = "royalblue") +
    labs(x = expression(theta), y = '事後の妥当性') +
    geom_area(data = filter(df, theta < .5), fill = 'royalblue')
area_1

この図で色付けされている部分の面積（現状では事後分布を標準化していないので、より正確には曲線下の面積全体に対する色付けされた部分の面積の比率）が、\(\theta < 0.5\) となる確率である。 この部分の面積を正確に求めるには、積分を計算する必要がある。 幸い、今回の問題のように簡単な例では、積分計算が可能である。

また、グリッド近似から直接面積を計算することもできる。 グリッドのうち、\(\theta < 0.5\) となる部分の確率をすべて足し合わせればよいので、次のコードで計算できる。

df %>% filter(theta < .5) %>% with(sum(posterior))
## [1] 0.1718746

この結果から、\(\theta < 0.5\) となる事後確率は 0.172であることがわかる。

しかし、いつもそう簡単に面積を求めることができるとは限らない。 そこで有効な方法が、サンプルを無作為に抽出し、そのサンプル内での特定の値の頻度を利用して計算する方法である。 シミュレートしたサンプルを使って、同様の計算を行ってみよう。

sum(post_sim < .5) / n_sim
## [1] 0.17422

この結果は、\(\theta < 0.5\) となる事後確率が 0.174 であることを示している。このように、シミュレートしたサンプルからでも、グリッド近似の結果から計算した結果と実質的にほぼ同じ結果が得られる。また、求め方は、サンプルを使う場合の方がやや簡単である。

次に、\(0.4 < \theta < 0.6\) となる確率を求めてみよう。 この確率を図示すると、次の図のようになる。

area_2 <- ggplot(df, aes(x = theta, y = posterior)) +
    geom_line(color = "royalblue") +
    labs(x = expression(theta), y = '事後の妥当性') +
    geom_area(data = filter(df, theta > .4 & theta < .6), fill = 'royalblue')
print(area_2)

この図の色付けされた部分の面積（正確には、全体の面積に対する色付けされた部分の面積の比率）が\(0.4 < \theta < 0.6\) となる確率である。シミュレートしたサンプルで計算してみよう。

sum(post_sim > .4 & post_sim < .6) / n_sim
## [1] 0.32973

この結果から、私たちは「Wの割合が 0.4から0.6 の間にある確率は 0.33 である」という（とりあえずの）結論を述べることができる。もちろん、他の区間、例えば \(0.7 < \theta < 0.9\) の確率や、\(0.2 < \theta 0.25\) の確率を提示することもできる。

分位数を求める

上の例では、ある一定の範囲を定め、\(\theta\) がその範囲に入る確率を計算した。 今度は、ある確率を定め、\(\theta\) の値がどの範囲にあれば、\(\theta\) がその範囲にある確率が指定した確率になるか考えよう。

まず、\(\theta\) の値が小さい方から考えて、確率（面積）が80% になる範囲を考えよう。 このような範囲は、シミュレートしたサンプルに対して quantile(x, probs) を使うことで簡単に求められる。 Quantile とは分位点のことである。今回は80パーセンタイルを知りたいので、probs = .8 を指定する。

(lower80 <- round(quantile(post_sim, probs = .8), 3))
##   80% 
## 0.762

このように、事後分布における\(\theta\) の80パーセンタイルは 0.762 であることがわかる。 つまり、\(0 < \theta <\) 0.762 となる確率は80パーセントである。

この面積を図示する。

area_l80 <- ggplot(df, aes(x = theta, y = posterior)) +
    geom_line(color = "royalblue") +
    geom_area(data = filter(df, theta < lower80), fill = 'royalblue') +
    labs(x = expression(theta), y = '事後の妥当性',
         title = '下位80パーセント')
print(area_l80)

次に、\(\theta\) 値が小さいもの10%と大きいもの10%を除いた中間の80%区間を求めよう。 この区間は、10パーセンタイルと90パーセンタイルの間なので、次のように求めることができる。

(mid80 <- round(quantile(post_sim, probs = c(.1, .9)), 3))
##   10%   90% 
## 0.446 0.813

よって、\(\theta \in\) 0.446, 0.813] となる確率が80パーセントである。

これを図示してみよう。

area_m80 <- ggplot(df, aes(x = theta, y = posterior)) +
    geom_line(color = "royalblue") +
    geom_area(data = filter(df, theta > mid80[1] & theta < mid80[2]), fill = 'royalblue') +
    labs(x = expression(theta), y = '事後の妥当性',
         title = '中間80パーセント')
print(area_m80)

これまでに求めた2つの区間は、いずれも\(\theta\) の80パーセントを表す区間であるが、区間の長さが異なる。 下側80パーセント区間より、中間80パーセント区間の方が短い。 これは、中間80パーセント区間の方が、確率密度がより高い範囲を捉えているためである。


最高事後確率密度区間 (HPDI)

上の例で、同じ確率（面積）を示す区間の取り方は複数あり、区間をどう取るかによって区間の長さが異なることが分かった。では、どのような区間を使って推論を行うべきだろうか。一つの考え方は、同じ確率を与える区間の中で、長さが最小になるものを選ぶというものである。区間の長さを最小にするには、確率密度が最大の部分から順番に面積を足し合わせていけばよい。このような方法で求められた区間のことを、最高事後確率密度区間 (highest posterior density interval: HPDI または HDI) と呼ぶ。

HPDI とその他の確率区間の違いをはっきりさせるために、地球儀問題で「3回トスして3回ともWが観察された」というデータを考える。

上で定義したグリッド近似の関数を使い、事後分布は次のように近似できる。

df_new <- globe_grid(w = 3, n = 3)
post_sim_new <- sample(df_new$theta, size = n_sim,
                       prob = df_new$posterior, replace = TRUE)

この事後分布の中間50パーセント区間は、分布の下側と上側からそれぞれ25パーセントずつ除外した部分なので、次のように求められる。

( mid50 <- quantile(post_sim_new, probs = c(.25, .75)) )
##       25%       75% 
## 0.7077077 0.9309309

この区間と50パーセントの面積は以下のように図示できる。

area_m50 <- ggplot(df_new, aes(x = theta, y = posterior)) +
    geom_line(color = "royalblue") +
    geom_area(data = filter(df_new, theta > mid50[1] & theta < mid50[2]),
              fill = 'royalblue') +
    labs(x = expression(theta), y = '事後の妥当性',
         title = '中間50パーセント')
print(area_m50)

このように、\(\theta\) が [0.708, 0.931] に入る事後確率が50パーセントであることがわかる。

しかし、この区間は、もっとも妥当性が高い \(\theta = 1\) 付近の値を除外してしまっている。 HPDI を使えば、より妥当性が高い部分を含んだ区間が求められる。

そこで、HPDIを計算する関数を以下のように作る。（現時点でこの関数の中身を理解できなくてもまったく問題ない。）

get_HPDI <- function(data, prob = .9) {
    # HPDI を計算する
    # 必要なパッケージ: ggplot2
    # 引数: data = theta と posterior を変数に持つデータフレーム
    #       prob = HPDIの大きさ、デフォルトは90パーセントHPDI
    # 返り値: リスト（HPDI, HPDIの図）
    
    # エラーチェック
    if (sum(c('theta', 'posterior') %in% names(data)) != 2) {
        stop(message = 'data must have columns "theta" and "posterior"')
    }
    if (prob < 0 | prob >1) stop(message = 'prob must be in [0, 1]')
    
    # HPDIの計算
    df <- data %>%
        arrange(desc(posterior)) %>%  # poserior が大きい順に並べ替え
        mutate(cum_sum = cumsum(posterior), # 累積確率を計算
               HPDI = ifelse(cum_sum <= prob, TRUE, FALSE)) # HDIに入るかどうか
    HPDI <- df %>%
        filter(HPDI == TRUE) %>%
        with(c(min(theta), max(theta)))
    
    # HPDIを視覚化
    plt <- ggplot(df, aes(x = theta, y = posterior)) +
        geom_line(color = 'firebrick') +
        geom_area(data = filter(df, HPDI == TRUE), fill = 'firebrick') +
        labs(x = expression(theta), y = '事後の妥当性',
             title = '事後分布とHPDI')
    return(list(HPDI = list(prob = prob, bounds = HPDI), plot = plt))
}

この関数を使って、50パーセントHPDI を求める。

hpdi50 <- get_HPDI(df_new, prob = .5)
hpdi50$HPDI
## $prob
## [1] 0.5
## 
## $bounds
## [1] 0.8418418 1.0000000

この区間を図示する。

hpdi50$plot

このように、\(\theta\) が[0.842, 1] に入る確率が50パーセントである。このHPDI は上で求めた中間50パーセント区間より短い。したがって、より精度の高い予測が行える。また、もっとも妥当性が高い\(\theta = 1\) 付近の値を区間に含んでいて、\(\theta\) の値として妥当性が高いのはどの範囲の値かを考えるの適している。

地球儀問題で「3回トスして3回ともWが観察された」というデータを分析すると、 「地球儀表面の水の割合が 84% 以上である確率は50%である」という推論ができる。

授業の内容に戻る