記事の目的
ポアソン分布と、その共役事前分布のガンマ分布を使用し、Rを使ってベイズ推定を行います
ある商品Aの平均購買個数を事後分布で推定します。
参考：ベイズ推論による機械学習入門
目次
0. モデルの説明
1. ライブラリ
2. 推定する分布
3. 事前分布
4. 事後分布
5. 予測分布
#0. モデルの説明
#1 ライブラリ
library(dplyr)
library(ggplot2)
set.seed(100)
#2. 推定する分布
商品Aの購買される個数は平均が15のポアソン分布に従います。しかし、僕たちはそれを知りません。
この、真の分布の平均の15を事後分布で推定します。
NULL %>% ggplot(aes(x=1:40)) + geom_bar(aes(y=dpois(1:40, 15)), stat="identity") +
labs(x="商品Aの購買個数", y="確率密度", title="推定する分布")
#3. 事前分布
事前分布として、ポアソン分布の共役事前分布であるガンマ分布を指定します。
データに適合しやすいように、a0, b0 は小さい値を初期値として設定します。
a0 <- 1
b0 <- 1
curve(dgamma(x, a0, b0), 0, 50 , xlab = "商品Aの平均購買個数(lambda)", ylab="確率密度", ylim=c(0,0.8))
#4. 事後分布
真の分布から50個のサンプルを取って事後分布を推定します。
事後分布は緑の曲線です。50付近が一番出やすいと推定できています。
#データ
N <- 50
X <- rpois(N, 15)
#事後分布
a <- sum(X)+a0
b <- N+b0
curve(dgamma(x, a, b), 0, 50 ,
xlab = "商品Aの平均購買個数(lambda)", ylab="確率密度", col="green", add=T)
#5. 予測分布
予測分布が、真の分布をうまく推定できていることが分かります。
Data <- data.frame(x=rep(1:40,2),
compare=rep(c("推定する分布", "予測分布"), each=40),
y = append(dpois(1:40, 15), dbinom(1:40, a, 1/(b+1))) )
Data %>% ggplot(aes(x=x, y=y, fill=compare)) +
geom_bar(stat="identity", alpha=0.7, position="identity") +
labs(x="商品Aの購買個数", y="確率密度", title="推定する分布と予測分布の比較") +
scale_fill_manual(values = c("black", "blue"))
