ベルヌーイ分布 : ベータ分布 (ベイズ推定)

記事の目的

ベルヌーイ分布と、その共役事前分布のベータ分布を使用し、Rを使ってベイズ推定を行います。
ある商品Aを購買する確率を推定します。
参考：ベイズ推論による機械学習入門

目次

0. モデルの説明
1. ライブラリ
2. 推定する分布
3. 事前分布
4. 事後分布
5. 予測分布

0. モデルの説明

1. ライブラリ

library(dplyr)
library(ggplot2)
set.seed(100)

2. 推定する分布

ある顧客が商品Aを購入する真の確率は0.8ですが、僕たちはそれを知りません。
この、真の確率0.8を事後分布で推定します。

NULL %>% ggplot(aes(x = c("1", "0"), y = c(0.8, 0.2))) + 
  geom_bar(stat = "identity") + ylim(0,1) + 
  labs(x="購買する(1) / 購買しない(0)", y="確率", title="推定する分布") +
  scale_x_discrete(limits=c("1", "0")) 

3. 事前分布

事前分布として、ベルヌーイ分布の共役事前分布であるベータ分布を指定します。
この事前分布は、商品Aを購入する確率は0.5が一番ありそうだが、あまり自信がないことを表しています。

a0 <- 2
b0 <- 2
curve(dbeta(x, a0, b0), 0,1, xlab="購買する確率", ylab="確率密度", ylim = c(0,10))

4. 事後分布

真の分布から100個サンプルを取って事後分布を推定します。
事後分布は緑の曲線です。0.8付近が一番出やすいと推定できています。

N <- 100
X <- rbinom(N, 1, 0.8)
a <- sum(X)
b <- N-a
curve(dbeta(x, a0+a, b0+b), 0,1, xlab="購買する確率", ylab="確率密度", add=T, col="green")

5. 予測分布

予測分布が、真の分布をうまく推定できていることが分かります。

Data <- data.frame(x=rep(c("1", "0"),2), y=c(0.8, 0.2, a/(a+b),b/(a+b)),
                   compare=rep(c("推定する分布", "予測分布"),each=2))
Data %>% ggplot(aes(x=x, y=y, fill=compare)) + 
  geom_bar(stat="identity", alpha=0.5, position="identity") +
  labs(x="購買する(1) / 購買しない(0)", y="確率密度", title="推定する分布と予測分布の比較") +
  scale_fill_manual(values = c("black", "blue")) +
  scale_x_discrete(limits=c("1", "0"))

◯SNS
・youtube
https://youtube.com/channel/UCFDyXEywtNhdtwqC3GAkYuA

・Twitter
https://twitter.com/Dken_ta