mcmcを使ってベイズ推定を自作でやってみる

ベイズ推定はライブラリやstanを使えばお手軽に実行可能
お手軽すぎて中身が本当に理解できているか不安
なので自作してみた

まずは準備

# 神ライブラリ
pacman::p_load(tidyverse)

# データの格納先
mu       <- NULL
lp       <- NULL

# パラメタ更新のステップサイズ
step     <- 0.01

データの用意

平均値0.5、分散1の正規分布に従う変数xを生成
今回はこのデータの平均値を推定してみる
推定値が0.5に近いければ実験成功

raw_data <- tibble(x = rnorm(n = 10000, mean = 0.5, sd = 1))

1つ目のiteration

推定するパラメタに初期値を与える
このときの対数尤度を求める

mu[1] <- runif(n = 1, min = 0, max = 1)

result  <- raw_data %>% 
  mutate(p = dnorm(x = x, mean = mu[1], sd = 1)) 

lp[1] <- sum(log(result$p))

2つ目以降のiteration

パラメタをステップサイズだけ更新
1. 更新の方向はランダムに選ぶ
このときの対数尤度を求める
更新前と更新後の対数尤度を比較して
1. {更新後 > 更新前} ならば、更新後のパラメタを採用して、1に戻る
2. {更新前 > 更新後} ならば
  1. {更新前の対数尤度 / 更新後の対数尤度}の比を求める
  2. その比の確率だけパラメタを更新するor更新しないを決めて、1に戻る

mu[2] <- mu[1] + step * if_else(runif(n = 1) >= 0.5, 1, -1)

for (i in 2:10000) {
  result  <- raw_data %>% 
    mutate(p = dnorm(x = x, mean = mu[i], sd = 1)) 
  
  lp[i] <- sum(log(result$p))

  if(lp[i] >= lp[(i - 1)]) {
    mu[(i + 1)] <- mu[i] + step * if_else(runif(n = 1) >= 0.5, 1, -1)
  } else {
    rate <- abs(lp[(i - 1)]) / abs(lp[i])
    mu_tmp <- sample(x = c(mu[(i - 1)], mu[i]), size = 1, prob = c(rate, 1 - rate))
    mu[(i + 1)] <- mu_tmp + step * if_else(runif(n = 1) >= 0.5, 1, -1)
  }
}

結果を見てみる

result <- tibble(id = 1:10000, mu = mu[1:10000], lp = lp[1:10000])

p1 <- ggplot(data = result) +
  geom_density(aes(x = mu), bw = 0.01) +
  geom_vline(aes(xintercept = mean(raw_data$x, na.rm = T)), col = 'red') +
  coord_cartesian(xlim = c(0, 1))

p2 <- ggplot(data = result) +
  geom_line(aes(x = id, y = mu))

patchwork::wrap_plots(p1, p2)

スクリプト

# 
pacman::p_load(tidyverse)

# 
raw_data <- tibble(x = rnorm(n = 10000, mean = 0.5, sd = 1))
mu       <- NULL
lp       <- NULL
step     <- 0.01

# 
mu[1] <- runif(n = 1, min = 0, max = 1)

result  <- raw_data %>% 
  mutate(p = dnorm(x = x, mean = mu[1], sd = 1)) 

lp[1] <- sum(log(result$p))

# 
mu[2] <- mu[1] + step * if_else(runif(n = 1) >= 0.5, 1, -1)

for (i in 2:10000) {
  
  result  <- raw_data %>% 
    mutate(p = dnorm(x = x, mean = mu[i], sd = 1)) 
  
  lp[i] <- sum(log(result$p))
  
  # 
  if(lp[i] >= lp[(i - 1)]) {
    mu[(i + 1)] <- mu[i] + step * if_else(runif(n = 1) >= 0.5, 1, -1)
  } else {
    rate <- abs(lp[(i - 1)]) / abs(lp[i])
    mu_tmp <- sample(x = c(mu[(i - 1)], mu[i]), size = 1, prob = c(rate, 1 - rate))
    mu[(i + 1)] <- mu_tmp + step * if_else(runif(n = 1) >= 0.5, 1, -1)
  }
}

result <- tibble(id = 1:10000, mu = mu[1:10000], lp = lp[1:10000])

p1 <- ggplot(data = result) +
  geom_density(aes(x = mu), bw = 0.01) +
  geom_vline(aes(xintercept = mean(raw_data$x, na.rm = T)), col = 'red') +
  coord_cartesian(xlim = c(0, 1))

p2 <- ggplot(data = result) +
  geom_line(aes(x = id, y = mu))

patchwork::wrap_plots(p1, p2)