準備

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

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