神戸ベイズ統計学研究会

Kruschke (2014), Ch.6

準備

必要なパッケージを読み込む。必要ならインストールする。

my_repos <- 'http://cran.ism.ac.jp/'
if (!require(ggplot2)) {
    install.packages('ggplot2', dependencies = TRUE, repos = my_repos)
    library('ggplot2')
}

## Loading required package: ggplot2

if (!require(dplyr)) {
    install.packages('dplyr', dependencies = TRUE, repos = my_repos)
    library('dplyr')
}

## Loading required package: 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

ベータ分布 (Beta Distribution)

ベータ分布は、尤度がベルヌーイ分布または二項分布のとき、共役事前分布 (conjugate prior) として利用される。 ベータ分布は2つの形状母数 (shape parameters) \(\alpha\) と \(\beta\) をもち、その確率密度関数は、 \[p(x | \alpha, \beta) = \frac{1}{B(\alpha, \beta)} x^{\alpha - 1} (1 - x)^{\beta - 1}\] である。ただし、 \[B(\alpha, \beta) = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha + \beta)} = \int_0^1 t^{\alpha - 1} (1 - t)^{\beta -1} dt\] である。また、確率変数 \(x\) は \(x\in[0, 1]\) で、\(\alpha\), \(\beta\) はともに正の実数である。 \(B(.,.)\) はベータ関数、\(\Gamma(.)\) はガンマ関数（階乗の一般化関数）である。

ベータ分布の平均値は、 \[\mu = \frac{\alpha}{\alpha + \beta}\] である。また、\(\alpha > 1\)、\(\beta > 1\) のとき、最頻値は \[\omega = \frac{\alpha - 1}{\alpha + \beta -2}\] となる。分散は、 \[\sigma^2 = \frac{\alpha \beta}{(\alpha + \beta)^2 (\alpha + \beta + 1)}\] である。

ベータ分布をよく理解するために、確率密度を図示する。 パラメタの値を変えて繰り返し使うために関数を作る。 返り値には、ベータ分布の特徴を表すリストを使う。

pdf_beta <- function(alpha, beta, var_name = expression(theta), fill = FALSE) {
    # Display the pdf of a beta distribution
    # Required packages: ggplot2
    # Args: alpha = first shape parameter, a positive real
    #       beta = second shape parameter, a positive real
    #       var_name = label for the random variable, default to theta
    #       fill = logical: set it TRUE to color the area under the curve
    # Return: list of parameters and stats
    if (alpha <= 0) stop(message = 'alpha must be a positive real number')
    if (beta <= 0) stop(message = 'beta must be a positive real number')
    theta <- seq(0, 1, length.out = 1000)
    density <- dbeta(theta, shape1 = alpha, shape2 = beta)
    df <- data.frame(theta = theta, density = density)
    plt <- ggplot(df, aes(theta, density)) +
        geom_line(color = 'darkblue') +
        labs(x = var_name, title = paste0('PDF of Beta(', alpha, ', ', beta, ')'))
    if (fill) plt <- plt + geom_area(fill = "darkblue")
    print(plt)
    mean <- alpha / (alpha + beta)
    mode <- (alpha - 1) / (alpha + beta - 2)
    sd <- sqrt( alpha*beta / ( (alpha + beta)^2 * (alpha + beta + 1) ) )
    lst <- list(param = list(shpae1 = alpha, shape2 = beta),
                stats = list(mean = mean, mode = mode, sd = sd),
                plot = plt)
    return(lst)
}

試しに使ってみる。

beta_test <- pdf_beta(alpha = 1, beta = 1)

beta_test <- pdf_beta(3, 3, var_name = 'x', fill = TRUE)

beta_test <- pdf_beta(12, 4, var_name = expression(gamma))

beta_test$stats

## $mean
## [1] 0.75
## 
## $mode
## [1] 0.7857143
## 
## $sd
## [1] 0.105021

beta_test <- pdf_beta(10, 20, fill = TRUE)

beta_test$stats$mode

## [1] 0.3214286

ベータ分布を事前分布として利用するとき、（超事前分布をおかないなら）2つの形状母数の値を決める必要がある。

例えば、特定の平均値 \(\mu \in (0, 1)\) に対応する形状母数は、 \[\alpha = \mu \kappa,\] \[\beta = (1 - \mu) \kappa\] から決められる。ただし、\(\kappa = \alpha + \beta\) で、この値は事前分布に観測値いくつ分の重みを与えるかを表す。

Kruschke (2014) がこの計算を実行する関数を用意してくれているが、練習のために自前で関数を作る。

beta_shapes_mean <- function(mean, kappa) {
    # Find the shape parameters of the beta that gives us a specific
    # mean and kappa
    # Args: mean in (0, 1)
    #       kappa > 0, the weight of the prior
    # Return: list of shape parameters
    if (mean <= 0 | mean >= 1) stop(message = 'mean must be between (0, 1)')
    if (kappa <=0) stop(message = 'kappa must be a positive real number')
    shape1 <- mean * kappa
    shape2 <- kappa - shape1
    return(list(shape1 = shape1, shape2 = shape2))
}

Kruscheke (2014:131) と同じ数値例で動作確認。

beta_shapes_mean(mean = .25, kappa = 4)    

## $shape1
## [1] 1
## 
## $shape2
## [1] 3

平均の代わりに、最頻値から形状母数を決めることもできる。 その場合、 \[\alpha = \omega (\kappa - 2) + 1,\] \[\beta = (1 - \omega)(\kappa - 2) + 1\] となる。ただし、ベータ分布の最頻値が定義されるのは\(\alpha >1\) かつ \(\beta\)>1 のときだけなので、\(\kappa >2\) である必要がある。これも自前の関数を作る。

beta_shapes_mode <- function(mode, kappa) {
    # Find the shape parameters of the beta that gives us a specific
    # mode and kappa
    # Args: mode in (0, 1)
    #       kappa > 2, the weight of the prior
    # Return: list of shape parameters
    if (mode <= 0 | mode >= 1) stop(message = 'mode must be between (0, 1)')
    if (kappa < 2) stop(message = 'kappa must be a real number larger than 2')
    shape1 <- mode * (kappa - 2) + 1
    shape2 <- kappa - shape1
    return(list(shape1 = shape1, shape2 = shape2))
}
beta_shapes_mode(mode = .25, kappa = 4)

## $shape1
## [1] 1.5
## 
## $shape2
## [1] 2.5

平均値と事前分布の重み \(\kappa\) の組み合わせの代わりに、平均値と標準偏差（分散）から形状母数を特定することもできる。 \[\alpha = \mu \left( \frac{\mu (1 - \mu)}{\sigma^2} - 1 \right),\] \[\beta = (1 - \mu) \left( \frac{\mu (1 - \mu)}{\sigma^2} - 1 \right)\] となる。 \[\sigma^2=\frac{\alpha \beta}{(\alpha+\beta)^2 (\alpha + \beta + 1)} = \frac{\mu(1 - \mu)}{\alpha + \beta + 1}\] なので、\(\sigma < 0.5\) である。実践的には、\(\alpha=\beta=1\) のときの標準偏差である \(0.289\) より小さい \(\sigma\) を使うことが多い。これについても関数を作る。

beta_shapes_mean_sd <- function(mean, sd) {
    # Find the shape parameters of the beta that gives us a specific
    # mean and a sd
    # Args: mean in (0, 1)
    #       sd in (0, 0.5)
    # Return: list of shape parameters
    if (mean <= 0 | mean >= 1) stop(message = 'mean must be between (0, 1)')
    if (sd <= 0 | sd >= .5) stop(message = 'sd must be between (0, 0.5)')
    shape1 <- mean * (mean * (1 - mean) / sd^2 - 1)
    shape2 <- (1 - mean) * (mean * (1 - mean) / sd^2 - 1)
    return(list(shape1 = shape1, shape2 = shape2))
}
beta_shapes_mean_sd(mean = .5, sd = .1)

## $shape1
## [1] 12
## 
## $shape2
## [1] 12

尤度

Kruschke (2014:Ch.6) では、尤度はベルヌーイ分布の積で表されている（二項分布で表現しても結果は同じだが）。 この尤度を図示する関数を作る。

データはすべて0または1でなければいけないので、まずはデータの不正をチェックする関数を作る。

check_bern <- function(data) {
    # Check if the data vector is the outcome of Bernoulli experiments
    # Arg: data = data vector
    # Return: TRUE if data contains only 0s and 1s, FALSE otherwise
    test <- data == 0 | data == 1
    return(length(test) == sum(test))
}

関数のテスト。

check_bern(c(0, 1, 0, 1, 1, 0))   # 正しいデータ

## [1] TRUE

check_bern(c(0, 1, 0, 1, .8, 0))  # 小数があるとき

## [1] FALSE

check_bern(c(0, 1, 0, 1, 1, 2))   # 0, 1 以外の整数があるとき

## [1] FALSE

これを利用して尤度関数を作る。

bern_lik <- function(data, var_name = expression(theta), fill = FALSE) {
    # Calculate the Bernoulli likelihood and display it
    # Required: ggplot2
    # Arg: data = vector of 0s and 1s
    #      var_name = label for the random variable, default to theta
    #      fill = logical: set it TRUE to color the area under the curve
    # Return: list of data frame containing some param values and likelihoods
    #                 and plt
    if (!check_bern(data)) stop(message = 'data values must be 0 or 1')
    z <- sum(data)
    n <- length(data)
    theta <- seq(0, 1, length.out = 1000)
    lik <- dbinom(z, size = n, prob = theta)
    df <- data.frame(param = theta, likelihood = lik)
    plt <- ggplot(df, aes(param, likelihood)) +
        geom_line(color = 'skyblue') +
        labs(x = var_name, title = 'Bernoulli Likelihood')
    if (fill) plt <- plt + geom_area(fill = 'skyblue')
    print(plt)
    return(list(data_frame = df, plot = plt))
}

テスト。

lik_test <- bern_lik(data = c(0))

lik_test <- bern_lik(data = c(1), fill = TRUE)

lik_test <- bern_lik(data = c(0, 1, 0, 1, 1), fill = TRUE)

事後分布

ベルヌーイ分布の尤度とベータ分布の事前分布から事後分布を求め、それを図示する関数を作る。 ついでに、HDI (highest density interval) を計算し、それを返り値にする。

post_bern_beta <- function(prior, data, prob_HDI = .90, fill = FALSE,
                           var_name = expression(theta), only_post = FALSE) {
    # Calculate posterior from Bernoulli likelihood and beta prior
    # Required: ggplot2 and dplyr
    # Args: prior = (1) 2-component vector of shape parameters of beta dist or
    #               (2) output of the last use of this function
    #       data = vector of 0s and 1s
    #       prob_HDI = prob for HDI, default to 90%
    #       fill = logical: TRUE to color the area under the curve
    #       var_name = label for the random variable, default to theta
    #       only_post = logical: TRUE if displays posterior only
    # Return: HDI
    if (prob_HDI <=0 | prob_HDI >=1) stop(message = 'prob_HDI must be in (0,1)')
    if (!check_bern(data)) stop(message = 'data values must be 0 or 1')
    if (class(prior) == "list"){
        if (is.null(prior$beta_shapes)) stop(message = 'invalid prior is used')
        prior <- prior$beta_shapes
    }
    if (length(prior) != 2) stop(message = '2-vector is required for prior')
    if (prior[1] <= 0 | prior[2] <= 0) {
        stop(message = 'each prior shape parameter must be a positive real number')
    }
    delta <- .0001
    theta <- seq(0, 1, by = delta)
    z <- sum(data)
    n <- length(data)
    alpha <- z + prior[1]
    beta <- n - z + prior[2]
    post <- dbeta(theta, shape1 = alpha, shape2 = beta)
    
    # Calculate HDI
    df <- data.frame(param = theta, post = post)    
    df <- df %>%
        arrange(desc(post)) %>%
        mutate(cum_prob = delta * cumsum(post),
               HDI = ifelse(cum_prob <= prob_HDI, TRUE, FALSE)) %>%
        arrange(param)
    HDI_interval <- df %>%
        filter(HDI == TRUE) %>%
        with(range(param))
    HDI <- list(percent = paste0(100*prob_HDI, "%"), interval = HDI_interval)
    
    # Visualizing posterior
    if (!only_post) {
        prior <- dbeta(theta, prior[1], prior[2])
        lik <- dbinom(z, size = n, prob = theta)
        df_large <- data.frame(
            param = rep(theta, 3),
            y = c(prior, lik, post),
            type = factor(rep(c(1, 2, 3), rep(length(theta), 3)),
                               labels = c('prior', 'likelihood', 'posterior')),
            HDI = c(rep(NA, 2*length(theta)), df$HDI)
        )
        plt <- ggplot(df_large, aes(param, y)) +
            geom_line(color = 'royalblue') +
            facet_grid(type ~ ., scales = 'free_y') +
            labs(x = var_name, y = "density / likelihood",
                 title = 'Bayesian Updating')
        if (fill) plt <- plt + geom_area(fill = 'royalblue')
    } else {
        plt <- ggplot(df) +
            geom_line(aes(param, post), color = "firebrick") +
            labs(x = var_name, y = "posterior density",
                 title = 'Posterior Distribution')
        if (fill) {
            df_HDI <- df %>%
                filter(HDI == TRUE) %>%
                mutate(zero = rep(0, n()))
            plt <- plt + geom_ribbon(data = df_HDI,
                                     aes(x = param, ymin = zero, ymax = post),
                                     fill = 'firebrick')
        }
    }
    print(plt)
    return(list(HDI = HDI, beta_shapes = c(alpha, beta), plot = plt))
}

試しに使ってみる（fill = TRUE にすると少し時間がかかる（かもしれない）ので注意）。

post_test <- post_bern_beta(prior = c(3, 4), data = c(1))

post_test <- post_bern_beta(prior = c(3, 4), data = c(1), fill = TRUE)

post_test$HDI

## $percent
## [1] "90%"
## 
## $interval
## [1] 0.2254 0.7746

post_test$beta_shapes

## [1] 4 4

fill = TRUE かつ post_only にすると、HDI が表示される。

# 90% HDI を表示
post_test <- post_bern_beta(prior = c(3, 4), data = c(1),
                            fill = TRUE, only_post = TRUE)

# 95% HDI を表示
post_test <- post_bern_beta(prior = c(3, 4), data = c(1), prob_HDI = .95,
                            fill = TRUE, only_post = TRUE)

練習問題

せっかくなので、自前の関数で解いていく。

Exercise 6.1

(A)

post_1 <-post_bern_beta(prior = c(4, 4), data = c(1), fill = TRUE)

(B)

post_2 <- post_bern_beta(prior = post_1, data = c(1), fill = TRUE)

(C)

post_3 <- post_bern_beta(post_2, data = c(0), fill = TRUE)

(D)

観測値1つ1つを順番に使ってアップデートを繰り返す。

post <- post_bern_beta(c(4, 4), data = c(0))

post <- post_bern_beta(post, data = c(1))

post <- post_bern_beta(post, data = c(1))

0を観察すると事後分布の山が左に移動し、1を観察すると右に移動することがわかる。

次に、データをまとめてアップデートする。

post <- post_bern_beta(c(4, 4), data = c(0, 1, 1))

このように、1つずつ更新しても、まとめて更新しても、最終的な結果（事後分布）はまったく同じものになる。

Exercise 6.2

(A)

作った関数に、問題の数字を与えれば解ける。

候補者A の支持率（予測得票率）を \(\theta\)、候補者B の支持率を \(1 - \theta\) とする。 また、事前分布が一様分布ということは、ベータ分布の形状母数を c(1, 1) にすればよい。 無作為に抽出された有権者100人のうち58人がAを、42人がBを好んだということなので、0が42個、1が58個のデータベクトルを 与える。

pred_elec <- post_bern_beta(prior = c(1, 1), data = rep(c(0, 1), c(42, 58)),
                            prob_HDI = .95, fill = TRUE, only_post = TRUE)

pred_elec$HDI

## $percent
## [1] "95%"
## 
## $interval
## [1] 0.4830 0.6731

95パーセントHDIは、\([0.48, 0.67]\) なので、どちらかといえば候補者A の方が有利であると言えそうだが、この結果だけで判断するのは難しい。

(B)

先ほどの結果を、新たな結果でアップデートする。新たに抽出された100人のうち、候補者A を好んだのは57人なので、次のようにアップデートする。

pred_elec <- post_bern_beta(prior = pred_elec, data = rep(c(0, 1), c(43, 57)),
                            prob_HDI = .95, fill = TRUE, only_post = TRUE)

pred_elec$HDI

## $percent
## [1] "95%"
## 
## $interval
## [1] 0.5062 0.6419

95パーセントHDIは、\([0.51, 0.64]\) なので、候補者Aが勝つ確率は95%以上である。

Exercise 6.3

この調査における結果変数は選ばれるキーであり、選択肢はFキーとJキーの 2つである。被験者 \(i\) がFキーを選ぶことを \(y_i =1\)、Jキーを選ぶことを \(y_i = 0\) とし、 Fが選ばれる確率を\(\theta \in [0, 1]\) とすると、 この問題の尤度はベルヌーイ分布でモデル化できる。 \[y_i | \theta \sim \mbox{Bern}(\theta),\] \[\Pr(y_i | \theta) = L_i(\theta | y_i) = \theta^y_i (1 - \theta)^{1-y_i}.\]

各被験者の選択が独立で、Fを好むバイアスは共通のもの（つまり、\(\theta\) は被験者間で同じ） とすると、全体の尤度は、 \[L(\theta | y) = \prod L_i(\theta | y_i)\] となる。したがって、\(n\) 人の被験者のうち \(z = \sum y_i\)人がFを選ぶときの尤度は、 \[L(\theta | y) = \theta^z (1 - \theta)^{n - z}\] である。

推定する母数は\(\theta\) であり、この\(\theta\) の事前確率を一様分布で表す。すなわち、 \[\theta \sim \mbox{Unif}(0, 1) = \mbox{Beta}(1, 1),\] \[p(\theta) = 1\] とする。

第1調査：2つのトレーニングに共通の単語 “radio” に対する反応

この調査の結果、\(n=50\)人のうち、\(z=40\)人がFを選んだので、事後分布は \[p(\theta | y) \propto \theta^{40} (1 - \theta)^{10}\] となる。よって、 \[\theta | y \sim \mbox{Beta}(41, 11)\] である。

これを図示する。

probe_1 <- post_bern_beta(prior = c(1, 1), data = rep(c(0, 1), c(10, 40)),
                          prob_HDI = .95, fill = TRUE, only_post = TRUE)

probe_1$HDI

## $percent
## [1] "95%"
## 
## $interval
## [1] 0.6773 0.8934

95パーセントHDIは\([0.68, 0.90]\)であり、0.5付近の値をとる確率は非常に小さいので、 F寄りのバイアスが存在していると言える。 すなわち、2つのトレーニングに共通する単語に対する反応では、最初に行ったトレーニングの影響の方が大きい。

第2調査：各トレーニングに固有の “ocean” と “moutain” を同時に示されたときの反応

この調査の結果、\(n=50\)人のうち、\(z=15\)人がFを選んだので、事後分布は \[p(\theta | y) \propto \theta^{15} (1 - \theta)^{35}\] となる。よって、 \[\theta | y \sim \mbox{Beta}(16, 36)\] である。

これを図示する。

probe_2 <- post_bern_beta(prior = c(1, 1), data = rep(c(0, 1), c(35, 15)),
                          prob_HDI = .95, fill = TRUE, only_post = TRUE)

probe_2$HDI

## $percent
## [1] "95%"
## 
## $interval
## [1] 0.1866 0.4329

95パーセントHDIは\([0.19, 0.43]\)であり、0.5付近の値をとる確率は非常に小さいので、 J寄りのバイアスが存在していると言える。 すなわち、それぞれのトレーニングに固有の単語を同時に見せられたときの反応については、後から行ったトレーニングの影響の方が大きい。

Exercise 6.4

どちらかの面が出易いマジックコインであることはほぼ間違いないが、表と裏のどちらに偏っているかはまったくわからないとする。このとき、コイン投げの結果は表が出る確率 \(\theta\) のベルヌーイ試行であると考え、\(\theta\) の事前確率をベータ分布で表す。

どちらに偏っているかわからないので、ベータ分布の平均を0.5 とする。また、どちらかに偏っているということは、分布の中央付近の値は取りにくいはずなので、標準偏差を0.499 とする（ベータ分布の標準偏差は最大でも0.5である）。

上で定義した関数を使って、この分布の形状母数を求める。

prior <- beta_shapes_mean_sd(mean = .5, sd = .499)
prior

## $shape1
## [1] 0.002006016
## 
## $shape2
## [1] 0.002006016

よって、事前分布の形状母数を2つとも0.02に設定し、その確率密度を図示する。

unusual_prior <- pdf_beta(alpha = .02, beta = .02, fill = TRUE)

このように、ベータ分布の形状母数を0に近い値にすると、\(\theta=0\) と \(\theta =1\) の近くだけ密度が高くなる事前分布が設定できる。

この事前確率を、「5回コインを投げたら4回表が出た」というデータ\(D\) で更新すると、 事後分布は、 \[\theta | D \sim \mbox{Beta}(4.02, 1.02)\] となる。

これを図示する。

post_coin <- post_bern_beta(c(.02, .02), data = c(0, rep(1, 4)), prob_HDI = .95,
                            fill = TRUE, only_post = TRUE)

post_coin$HDI

## $percent
## [1] "95%"
## 
## $interval
## [1] 0.4706 0.9999

95パーセントHDIは\([0.47, 0.9999]\)であり、最頻値は\(\theta = 0.9999\) なので、表が出易い細工が施されたコインである可能性が高い。

Exercise 6.5

コインを1回投げて表がでる確率を \(\theta\) とする。 また、尤度はベルヌーイ分布でモデル化する。

(A)

\(\theta = 0.5\) であるという強い信念を持っていて、この信念は観測値 2,000個分に相当すると仮定する。

prior <- beta_shapes_mean(mean = .5, kappa = 2000)
prior

## $shape1
## [1] 1000
## 
## $shape2
## [1] 1000

となるから、この信念は、\(\theta \sim \mbox{Beta}(1000, 1000)\) で表現できる。

この事前確率を、「10回コインを投げて9回表」というデータ\(D\) で更新すると、事後分布は \[\theta | D \sim \mbox{Beta}(1009, 1001)\] となる。

これを図示すると、

post_strong_belief <- post_bern_beta(c(1000, 1000), data = c(0, rep(1, 9)), 
                                     fill = TRUE, only_post = TRUE)

このように、非常に強い信念（1点に集中した事前確率）からスタートすると、10個程度の観測値の影響はほとんどない。 したがって、11回目のコイン投げで表が出る確率は、0.5 であると予測できる。

(B)

マジックショップ製のコインであり、表か裏のどちらかに偏っていると思われるの、 先ほどの Ex.6.4 で使ったのと同じ事前分布、\(\theta \sim \mbox{Beta}(0.02, 0.02)\) を利用する。

この事前確率を、「10回コインを投げて9回表」というデータ\(D\) で更新すると、事後分布は \[\theta | D \sim \mbox{Beta}(9.02, 1.02)\] となる。

これを図示すると、

post_magic <- post_bern_beta(c(.02, .02), data = c(0, rep(1, 9)), 
                             fill = TRUE, only_post = TRUE)

post_magic$HDI

## $percent
## [1] "90%"
## 
## $interval
## [1] 0.7717 0.9999

このように、事前分布に与えられるウェイト（ベータ分布の形状母数の値の合計）が小さいと、事後分布はデータに支配される。

この事後分布の平均値は、

9.02 / (9.02 + 1.02)

## [1] 0.8984064

なので、11回目に表がでる確率は、0.9であると予測する。

Go Back