ベイズ統計02-二項分布とベータ分布-
UPDATE: 2024-01-07 16:08:52.744396
はじめに
このノートは「ベイズ統計」に関する何らかの内容をまとめ、ベイズ統計への理解を深めていくために作成している。今回は「社会科学のためのベイズ統計モデリング」の第8章を写経していく。基本的には気になった部分を写経しながら、ところどころ自分用の補足をメモすることで、自分用の補足資料になることを目指す。私の解釈がおかしく、メモが誤っている場合があるので注意。
8.2.1 藤井七段の対戦成績
今回使用するデータは将棋棋士の藤井七段の100試合の勝敗データ。この100試合のデータから藤井七段の強さを推定することが今回の分析目的。単純集計だと85勝15敗で勝率85%となる。
library(dplyr)
library(rstan)
library(ggplot2)
options(max.print = 999999)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())
d <- read.table('https://raw.githubusercontent.com/HiroshiHamada/BMS/master/ch08/Fujii.txt')
table(d)## V1
## 0 1
## 15 85# ggplot(data = data,aes(x = game,y = win+lose,fill = win)) +
# geom_point(size = 2.5, shape = 21, stroke = 0.5) +
# scale_fill_gradient(low = "black",high = "white") +
# scale_x_continuous(breaks = c(1,seq(10,100,5))) +
# ylim(0.9,1.1) +
# theme( panel.background = element_blank(),
# panel.grid = element_blank(),
# axis.title.x = element_blank(),
# axis.text = element_blank(),
# axis.title.y = element_blank(),
# axis.text.y = element_blank(),
# axis.ticks.y = element_blank(),
# legend.position = "none") 8.2.2 ベルヌーイ・モデル
ここでは、勝ちを1、負けを0とする確率変数\(Y_{i}\)とする。\(Y_{i}\)は独立であり、\(q\)をパラメタとするベルヌーイ分布に従うと仮定する。事前分布は\(Beta(\alpha=1,\beta=1)\)のベータ分布に従うと仮定する。
\[ \begin{eqnarray} Y_{i} &\sim& Bernoulli(q), \ \ i = 1...n\\ q &\sim& Beta(\alpha, \beta) \end{eqnarray} \]
ここで使用するベルヌーイ・モデルについては、第3章でわかりやすく解説されている。内容をメモしたものを貼り付けておく。
\(n+1\)試合目の予測分布は下記のパラメタを持つベルヌーイ分布に従う。
\[ E[q] = \frac{\alpha+\sum y_{i}}{\alpha+\beta+n} \]
8.2.3 藤井七段のデータ分析
勝利確率\(q\)の事後分布と予測分布を1試合終るごとにデータから推定する。データを得る前の状態は不明なので、\(q\)の事前分布は\(Beta(\alpha=1,\beta=1)\)のベータ分布に従うと仮定する。
data <- tibble(
game = 1:100,
win = d$V1,
lose = 1-win
) %>%
mutate(
cumwin = cumsum(win),
cumlose = cumsum(lose),
mle = cumwin/game
)
print(data, n = 100)## # A tibble: 100 × 6
## game win lose cumwin cumlose mle
## <int> <int> <dbl> <int> <dbl> <dbl>
## 1 1 1 0 1 0 1
## 2 2 1 0 2 0 1
## 3 3 1 0 3 0 1
## 4 4 1 0 4 0 1
## 5 5 1 0 5 0 1
## 6 6 1 0 6 0 1
## 7 7 1 0 7 0 1
## 8 8 1 0 8 0 1
## 9 9 1 0 9 0 1
## 10 10 1 0 10 0 1
## 11 11 1 0 11 0 1
## 12 12 1 0 12 0 1
## 13 13 1 0 13 0 1
## 14 14 1 0 14 0 1
## 15 15 1 0 15 0 1
## 16 16 1 0 16 0 1
## 17 17 1 0 17 0 1
## 18 18 1 0 18 0 1
## 19 19 1 0 19 0 1
## 20 20 1 0 20 0 1
## 21 21 1 0 21 0 1
## 22 22 1 0 22 0 1
## 23 23 1 0 23 0 1
## 24 24 1 0 24 0 1
## 25 25 1 0 25 0 1
## 26 26 1 0 26 0 1
## 27 27 1 0 27 0 1
## 28 28 1 0 28 0 1
## 29 29 1 0 29 0 1
## 30 30 0 1 29 1 0.967
## 31 31 1 0 30 1 0.968
## 32 32 1 0 31 1 0.969
## 33 33 0 1 31 2 0.939
## 34 34 1 0 32 2 0.941
## 35 35 1 0 33 2 0.943
## 36 36 1 0 34 2 0.944
## 37 37 0 1 34 3 0.919
## 38 38 1 0 35 3 0.921
## 39 39 1 0 36 3 0.923
## 40 40 1 0 37 3 0.925
## 41 41 1 0 38 3 0.927
## 42 42 0 1 38 4 0.905
## 43 43 0 1 38 5 0.884
## 44 44 1 0 39 5 0.886
## 45 45 0 1 39 6 0.867
## 46 46 1 0 40 6 0.870
## 47 47 1 0 41 6 0.872
## 48 48 1 0 42 6 0.875
## 49 49 0 1 42 7 0.857
## 50 50 1 0 43 7 0.86
## 51 51 1 0 44 7 0.863
## 52 52 1 0 45 7 0.865
## 53 53 1 0 46 7 0.868
## 54 54 1 0 47 7 0.870
## 55 55 1 0 48 7 0.873
## 56 56 1 0 49 7 0.875
## 57 57 0 1 49 8 0.860
## 58 58 1 0 50 8 0.862
## 59 59 1 0 51 8 0.864
## 60 60 0 1 51 9 0.85
## 61 61 1 0 52 9 0.852
## 62 62 1 0 53 9 0.855
## 63 63 1 0 54 9 0.857
## 64 64 0 1 54 10 0.844
## 65 65 1 0 55 10 0.846
## 66 66 0 1 55 11 0.833
## 67 67 1 0 56 11 0.836
## 68 68 1 0 57 11 0.838
## 69 69 1 0 58 11 0.841
## 70 70 1 0 59 11 0.843
## 71 71 1 0 60 11 0.845
## 72 72 1 0 61 11 0.847
## 73 73 1 0 62 11 0.849
## 74 74 1 0 63 11 0.851
## 75 75 1 0 64 11 0.853
## 76 76 1 0 65 11 0.855
## 77 77 1 0 66 11 0.857
## 78 78 1 0 67 11 0.859
## 79 79 1 0 68 11 0.861
## 80 80 1 0 69 11 0.862
## 81 81 1 0 70 11 0.864
## 82 82 1 0 71 11 0.866
## 83 83 0 1 71 12 0.855
## 84 84 1 0 72 12 0.857
## 85 85 1 0 73 12 0.859
## 86 86 1 0 74 12 0.860
## 87 87 1 0 75 12 0.862
## 88 88 1 0 76 12 0.864
## 89 89 1 0 77 12 0.865
## 90 90 1 0 78 12 0.867
## 91 91 0 1 78 13 0.857
## 92 92 1 0 79 13 0.859
## 93 93 1 0 80 13 0.860
## 94 94 0 1 80 14 0.851
## 95 95 1 0 81 14 0.853
## 96 96 0 1 81 15 0.844
## 97 97 1 0 82 15 0.845
## 98 98 1 0 83 15 0.847
## 99 99 1 0 84 15 0.848
## 100 100 1 0 85 15 0.851試合が終わるごとに事後分布と事後予測分布を計算することで、強さの推移を可視化する。1試合が終わるごとに計算する内容を書き下すと下記の通りである。
- 00試合目の結果より00勝00敗、予測分布による01試合目の予測勝利確率は\(\frac{1+0}{1+1+0} = 0.500\)となる
- 01試合目の結果より01勝00敗、予測分布による02試合目の予測勝利確率は\(\frac{1+1}{1+1+1} = 0.667\)となる
- 02試合目の結果より02勝00敗、予測分布による03試合目の予測勝利確率は\(\frac{1+2}{1+1+2} = 0.750\)となる
- 03試合目の結果より03勝00敗、予測分布による04試合目の予測勝利確率は\(\frac{1+3}{1+1+3} = 0.800\)となる
- 04試合目の結果より04勝00敗、予測分布による05試合目の予測勝利確率は\(\frac{1+3}{1+1+3} = 0.833\)となる
- 05試合目の結果より05勝00敗、予測分布による06試合目の予測勝利確率は\(\frac{1+4}{1+1+4} = 0.857\)となる
- …snip
- 29試合目の結果より29勝0敗、予測分布による30試合目の予測勝利確率は\(\frac{1+29}{1+1+29} = 0.968\)となる
- 30試合目の結果より29勝1敗、予測分布による31試合目の予測勝利確率は\(\frac{1+29}{1+1+30} = 0.938\)となる
- …snip
- 100試合目の結果より29勝1敗、予測分布による101試合目の予測勝利確率は\(\frac{1+85}{1+1+100} = 0.843\)となる
可視化するためのデータを用意する。事後分布の平均、95%ベイズ信頼区間を計算している。
## posterior
postdata <- tibble(
game = 0:100,
win = c(NA, data$win),
lose = c(NA, data$lose),
cumwin = c(NA,data$cumwin),
cumlose = c(NA,data$cumlose),
post_alpha = c(1,data$cumwin+1),
post_beta = c(1,data$cumlose+1),
mle = c(NA,data$mle),
post_mu = post_alpha / (post_alpha+post_beta),
post_med = qbeta(0.500, post_alpha, post_beta),
CI_lower = qbeta(0.025, post_alpha, post_beta),
CI_upper = qbeta(0.975, post_alpha, post_beta)
)
print(postdata, n = 100)## # A tibble: 101 × 12
## game win lose cumwin cumlose post_alpha post_beta mle post_mu
## <int> <int> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 NA NA NA NA 1 1 NA 0.5
## 2 1 1 0 1 0 2 1 1 0.667
## 3 2 1 0 2 0 3 1 1 0.75
## 4 3 1 0 3 0 4 1 1 0.8
## 5 4 1 0 4 0 5 1 1 0.833
## 6 5 1 0 5 0 6 1 1 0.857
## 7 6 1 0 6 0 7 1 1 0.875
## 8 7 1 0 7 0 8 1 1 0.889
## 9 8 1 0 8 0 9 1 1 0.9
## 10 9 1 0 9 0 10 1 1 0.909
## 11 10 1 0 10 0 11 1 1 0.917
## 12 11 1 0 11 0 12 1 1 0.923
## 13 12 1 0 12 0 13 1 1 0.929
## 14 13 1 0 13 0 14 1 1 0.933
## 15 14 1 0 14 0 15 1 1 0.938
## 16 15 1 0 15 0 16 1 1 0.941
## 17 16 1 0 16 0 17 1 1 0.944
## 18 17 1 0 17 0 18 1 1 0.947
## 19 18 1 0 18 0 19 1 1 0.95
## 20 19 1 0 19 0 20 1 1 0.952
## 21 20 1 0 20 0 21 1 1 0.955
## 22 21 1 0 21 0 22 1 1 0.957
## 23 22 1 0 22 0 23 1 1 0.958
## 24 23 1 0 23 0 24 1 1 0.96
## 25 24 1 0 24 0 25 1 1 0.962
## 26 25 1 0 25 0 26 1 1 0.963
## 27 26 1 0 26 0 27 1 1 0.964
## 28 27 1 0 27 0 28 1 1 0.966
## 29 28 1 0 28 0 29 1 1 0.967
## 30 29 1 0 29 0 30 1 1 0.968
## 31 30 0 1 29 1 30 2 0.967 0.938
## 32 31 1 0 30 1 31 2 0.968 0.939
## 33 32 1 0 31 1 32 2 0.969 0.941
## 34 33 0 1 31 2 32 3 0.939 0.914
## 35 34 1 0 32 2 33 3 0.941 0.917
## 36 35 1 0 33 2 34 3 0.943 0.919
## 37 36 1 0 34 2 35 3 0.944 0.921
## 38 37 0 1 34 3 35 4 0.919 0.897
## 39 38 1 0 35 3 36 4 0.921 0.9
## 40 39 1 0 36 3 37 4 0.923 0.902
## 41 40 1 0 37 3 38 4 0.925 0.905
## 42 41 1 0 38 3 39 4 0.927 0.907
## 43 42 0 1 38 4 39 5 0.905 0.886
## 44 43 0 1 38 5 39 6 0.884 0.867
## 45 44 1 0 39 5 40 6 0.886 0.870
## 46 45 0 1 39 6 40 7 0.867 0.851
## 47 46 1 0 40 6 41 7 0.870 0.854
## 48 47 1 0 41 6 42 7 0.872 0.857
## 49 48 1 0 42 6 43 7 0.875 0.86
## 50 49 0 1 42 7 43 8 0.857 0.843
## 51 50 1 0 43 7 44 8 0.86 0.846
## 52 51 1 0 44 7 45 8 0.863 0.849
## 53 52 1 0 45 7 46 8 0.865 0.852
## 54 53 1 0 46 7 47 8 0.868 0.855
## 55 54 1 0 47 7 48 8 0.870 0.857
## 56 55 1 0 48 7 49 8 0.873 0.860
## 57 56 1 0 49 7 50 8 0.875 0.862
## 58 57 0 1 49 8 50 9 0.860 0.847
## 59 58 1 0 50 8 51 9 0.862 0.85
## 60 59 1 0 51 8 52 9 0.864 0.852
## 61 60 0 1 51 9 52 10 0.85 0.839
## 62 61 1 0 52 9 53 10 0.852 0.841
## 63 62 1 0 53 9 54 10 0.855 0.844
## 64 63 1 0 54 9 55 10 0.857 0.846
## 65 64 0 1 54 10 55 11 0.844 0.833
## 66 65 1 0 55 10 56 11 0.846 0.836
## 67 66 0 1 55 11 56 12 0.833 0.824
## 68 67 1 0 56 11 57 12 0.836 0.826
## 69 68 1 0 57 11 58 12 0.838 0.829
## 70 69 1 0 58 11 59 12 0.841 0.831
## 71 70 1 0 59 11 60 12 0.843 0.833
## 72 71 1 0 60 11 61 12 0.845 0.836
## 73 72 1 0 61 11 62 12 0.847 0.838
## 74 73 1 0 62 11 63 12 0.849 0.84
## 75 74 1 0 63 11 64 12 0.851 0.842
## 76 75 1 0 64 11 65 12 0.853 0.844
## 77 76 1 0 65 11 66 12 0.855 0.846
## 78 77 1 0 66 11 67 12 0.857 0.848
## 79 78 1 0 67 11 68 12 0.859 0.85
## 80 79 1 0 68 11 69 12 0.861 0.852
## 81 80 1 0 69 11 70 12 0.862 0.854
## 82 81 1 0 70 11 71 12 0.864 0.855
## 83 82 1 0 71 11 72 12 0.866 0.857
## 84 83 0 1 71 12 72 13 0.855 0.847
## 85 84 1 0 72 12 73 13 0.857 0.849
## 86 85 1 0 73 12 74 13 0.859 0.851
## 87 86 1 0 74 12 75 13 0.860 0.852
## 88 87 1 0 75 12 76 13 0.862 0.854
## 89 88 1 0 76 12 77 13 0.864 0.856
## 90 89 1 0 77 12 78 13 0.865 0.857
## 91 90 1 0 78 12 79 13 0.867 0.859
## 92 91 0 1 78 13 79 14 0.857 0.849
## 93 92 1 0 79 13 80 14 0.859 0.851
## 94 93 1 0 80 13 81 14 0.860 0.853
## 95 94 0 1 80 14 81 15 0.851 0.844
## 96 95 1 0 81 14 82 15 0.853 0.845
## 97 96 0 1 81 15 82 16 0.844 0.837
## 98 97 1 0 82 15 83 16 0.845 0.838
## 99 98 1 0 83 15 84 16 0.847 0.84
## 100 99 1 0 84 15 85 16 0.848 0.842
## # ℹ 1 more row
## # ℹ 3 more variables: post_med <dbl>, CI_lower <dbl>, CI_upper <dbl>連勝している間は1と推定する最尤推定と比べ、事後分布は不確実性を考慮して、最初の方は信用区間も広く、事後平均も低い。多くのデータが観察できた後半の試合では、平均も安定し、信用区間も狭くなっている。
ggplot(data = postdata) +
geom_ribbon(aes(x = game,y = post_mu,
ymin = CI_lower, ymax = CI_upper), fill = 'grey80')+
geom_line(aes(x = game,y = mle,linetype = 'MLE')) +
geom_line(aes(x = game,y = post_mu, linetype= 'post. dist.')) +
labs(linetype = 'Estimation') +
ggtitle('Change in posterior mean of posterior distribution') +
theme_bw() 
事後分布の確率密度関数の変化をプロットする。最初は先程と同じく不確実性を考慮しているものの、最後の方は安定した分布になっている。
# x <- seq(0,1,0.001)
# plot(x, dbeta(x, shape1 = 86, shape2 = 16), type = 'l')
ggplot(data = tibble(q = c(0,1)), aes(x = q)) +
purrr::pmap(list(a = postdata$post_alpha, b = postdata$post_beta, co = postdata$game),
function(a, b, co)
stat_function(fun = dbeta, args = list(shape1 = a, shape2 = b), aes_q(color = co), size = 0.1)
) + labs(color = 'game') +
scale_colour_gradient(low = 'grey80', high = 'black') +
scale_x_continuous(breaks = seq(0, 1, 0.05)) +
ggtitle('Change in probability density function of posterior distribution') +
theme_bw()
8.2.5 先手後手で強さが変わるのか
一般的には将棋は先手が有利で、後手が不利とも言われる。先手後手の要因を入れてモデルを拡張する。\(i\)試合目の先手後手を区別する変数\(x_{i}\)を導入し、先手は\(x=1\)とする。
\[ \begin{eqnarray} Y_{i} &\sim& Bernoulli(x_{i}q_{1} + (1-x_{i})q_{0}), \ \ i = 1...n\\ q_{0} &\sim& Beta(\alpha_{0}, \beta_{0}) \\ q_{1} &\sim& Beta(\alpha_{1}, \beta_{1}) \\ \end{eqnarray} \]
Stanのモデルは下記の通り。
data {
int N;
int X[N];
int Z[N];
}
parameters {
real<lower=0,upper=1> q1;
real<lower=0,upper=1> q0;
}
model {
for (n in 1:N) {
X[n] ~ bernoulli(Z[n]*q1+(1-Z[n])*q0);
}
}このモデルを少し深掘りしておく。下記の通り、zが切り替わることで、x[n]が従う分布が異なるパラメタから生成される考えている。
model {
for (n in 1:N) {
X[n] ~ bernoulli(Z[n]*q1+(1-Z[n])*q0);
}
}
// z = 1のとき
// X[n] ~ bernoulli(1*q1); -> X[n] ~ bernoulli(q1);
// z = 0のとき
// X[n] ~ bernoulli(Z(1-0)*q0); -> X[n] ~ bernoulli(q0);データを用意して、
d <- read.table('https://raw.githubusercontent.com/HiroshiHamada/BMS/master/ch08/Fujii_first.txt', header = TRUE)
data <- list(N = nrow(d), X = d$win, Z = d$first)
data## $N
## [1] 100
##
## $X
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 1 0
## [38] 1 1 1 1 0 0 1 0 1 1 1 0 1 1 1 1 1 1 1 0 1 1 0 1 1 1 0 1 0 1 1 1 1 1 1 1 1
## [75] 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 0 1 0 1 1 1 1
##
## $Z
## [1] 0 1 0 1 1 1 0 0 1 1 0 1 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 1 1 0 1 0 0 1 0 1 0
## [38] 0 0 0 0 0 0 0 0 1 1 1 1 0 1 0 1 0 0 0 0 0 1 1 1 0 0 1 1 1 0 0 1 1 1 0 0 0
## [75] 1 1 1 0 1 1 1 1 0 1 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 1ここでは、stan_model()関数で最初にコンパイルしておいてから、
model <- stan_model('note_bayes02−001.stan')sampling()関数でサンプリングする。
fit <- sampling(object = model, data = data, seed = 1989)推定結果を確認する。後手と信用区間がすこし被るものの、先手の方が勝率が高くなっているため、藤井七段については先手の方が有利と言えそうである。
- 先手の場合、勝率は
q1=0.89[0.79-0.96] - 後手の場合、勝率は
q1=0.80[0.69-0.84]
print(fit)## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## q1 0.89 0.00 0.05 0.79 0.86 0.89 0.92 0.96 2584 1
## q0 0.80 0.00 0.05 0.69 0.77 0.80 0.84 0.89 3047 1
## lp__ -46.55 0.03 1.06 -49.32 -46.99 -46.21 -45.79 -45.52 1696 1
##
## Samples were drawn using NUTS(diag_e) at Sun Jan 7 16:09:18 2024.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).事後分布を可視化しておく。
stan_plot(
fit,
point_est = "mean",
ci_level = 0.95,
outer_level = 1.00,
show_density = TRUE,
fill_color = "grey") +
theme_bw()