# 6章章末問題
## 3 Model improvement
### Q
(a) Use the solution to the previous problem and your substantive knowledge to construct an improved model for airline fatalities.
(b) Fit the new model to the airline fatality data.
\(c\) Use your new model to forecast the airline fatalities in 1986. How does this differ from the forecasts from the previous models?
(d) Check the new model using the same posterior predictive checks as you used in the previous models. Does the new model fit better?
## 5
### Question
Hypothesis testing: discuss the statement, ‘Null hypotheses of no difference are usually known to be false before the data are collected; when they are, their rejection or acceptance simply reflects the size of the sample and the power of the test, and is not a contribution to science' (Savage, 1957, quoted in Kish, 1965). If you agree with this statement, what does this say about the model checking discussed in this chapter?
### Solution
サンプルサイズや検定力故に、仮設検定の意義が低い状況においては、その不確かさを事後分布の形で示すと良い。それにより、検定の棄却/採択の二者択一によるリスクを減らすことができる。
## 9
### Question
Model checking: check the assumed model fitted to the rat tumor data in Section 5.3. Define some test quantities that might be of scientific interest, and compare them to their posterior predictive distributions.
### Solution
元データは、70の実験($j \in \{ 1, \cdots 70 \}$)についてそれぞれ$\theta_j$を仮定し、それぞれ一定数($n_j$)のネズミの腫瘍の発生件数($y_j$)を調べたもの。
この問題は以下の順番で解いていくことができる。
1 $p(\alpha, \beta|\mathbf{y})$を導出
2 $p(\theta|\alpha, \beta)$を導出
3 70の実験についてそれぞれ$\theta_j$をサンプリングし、$y_j$をサンプリングし、$y_j/n_j$、試験量を導出
4 以上を200回繰り返し、実際のデータと比較することによって、$\alpha, \beta$の(そして、それを前提としたモデルの)妥当性を検証する。
1については、5章の議論により
\begin{align}
p(\alpha, \beta|\mathbf{y}) \propto p(\alpha, \beta) \prod_j \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\alpha)}\frac{\Gamma(\alpha+y_j)\Gamma(\beta+n_j-y_j)}{\Gamma(\alpha+\beta+n_j)}
\end{align}
である。このことは(解析的に計算することは困難だが)$p(\alpha,\beta)$の事前分布さえ設定すれば、$p(\alpha, \beta|\mathbf{y})$を導出することは、ガンマ関数をそれぞれ計算することによって容易にできることを意味する。
以下はRで。
```{r efwfwwf}
library(ggplot2)
theme_set(theme_minimal())
library(gridExtra)
library(tidyr)
library(latex2exp)
y <- c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,
2,1,5,2,5,3,2,7,7,3,3,2,9,10,4,4,4,4,4,4,4,10,4,4,4,5,11,12,
5,5,6,5,6,6,6,6,16,15,15,9,4)
n <- c(20,20,20,20,20,20,20,19,19,19,19,18,18,17,20,20,20,20,19,19,18,18,25,24,
23,20,20,20,20,20,20,10,49,19,46,27,17,49,47,20,20,13,48,50,20,20,20,20,
20,20,20,48,19,19,19,22,46,49,20,20,23,19,22,20,20,20,52,46,47,24,14)
yn_obs <- y/n
t_st <- data.frame(
t_rep = rep(NA, 200),
t_mean = rep(NA, 200)
)
#n番目を返す関数
n_order <- function(x, n){
sort(x )[n]
}
#まず、p(\alpha, \beta|y)を、\alpha, \betaごとに求める。
#https://avehtari.github.io/BDA_R_demos/demos_ch5/demo5_1.html 参照
A <- seq(0.5, 6, length.out = 100)
B <- seq(3, 33, length.out = 100)
cA <- rep(A, each = length(B))
cB <- rep(B, length(A))
lpfun <- function(a, b, y, n){
log(a+b)*(-5/2) + sum(lgamma(a+b)-lgamma(a)-lgamma(b)+lgamma(a+y)+lgamma(b+n-y)-lgamma(a+b+n))
}
lp <- mapply(lpfun, cA, cB, MoreArgs = list(y, n))
df_marg <- data.frame(x = cA, y = cB, p = exp(lp - max(lp)))
title1 <- TeX('The marginal of $\\alpha$ and $\\beta$')
ggplot(data = df_marg, aes(x = x, y = y)) +
geom_raster(aes(fill = p, alpha = p), interpolate = T) +
geom_contour(aes(z = p), colour = 'black', size = 0.2) +
coord_cartesian(xlim = c(1,5), ylim = c(4, 26)) +
labs(x = TeX('$\\alpha$'), y = TeX('$\\beta$'), title = title1) +
scale_fill_gradient(low = 'yellow', high = 'red', guide = F) +
scale_alpha(range = c(0, 1), guide = F)
#ここからthetaの標本を抽出
nsamp <- 100
samp_indices <- sample(length(df_marg$p), size = nsamp,
replace = T, prob = df_marg$p/sum(df_marg$p))
samp_A <- cA[samp_indices[1:nsamp]]
samp_B <- cB[samp_indices[1:nsamp]]
df_psamp <- mapply(function(a, b, x) dbeta(x, a, b),
samp_A, samp_B, MoreArgs = list(x = x)) %>%
as.data.frame() %>% cbind(x) %>% gather(ind, p, -x)
indtonum <- function(x) strtoi(substring(x,2))
title2 <- TeX('Beta($\\alpha,\\beta$) given posterior draws of $\\alpha$ and $\\beta$')
plot_psamp <- ggplot(data = subset(df_psamp, indtonum(ind) <= 20)) +
geom_line(aes(x = x, y = p, group = ind), color='forestgreen') +
labs(x = expression(theta), y = '', title = title2) +
scale_y_continuous(breaks = NULL)
df_psampmean <- spread(df_psamp, ind, p) %>% subset(select = -x) %>%
rowMeans() %>% data.frame(x = x, p = .)
title3 <- TeX('Population distribution (prior) for $\\theta_j$')
plot_psampmean <- ggplot(data = df_psampmean) +
geom_line(aes(x = x, y = p), color='forestgreen') +
labs(x = expression(theta), y = '', title = title3) +
scale_y_continuous(breaks = NULL)
grid.arrange(plot_psamp, plot_psampmean)
df_psampmean
#ここでは簡単に、平均的な分布から70個thetaを抽出し、それぞれについてy/nを求め、61番目と6番目に大きいy/nの差と、y/nをt統計量とし、それを200回繰り返す
#シミュレーション200回
for (i in 1:200) {
#実験毎のthetaをサンプリング
yn_rep <- rep(NA,70)
for (j in 1:70) {
theta_rep <- sample(df_psampmean$x, size = 1, replace = T, prob =df_psampmean$p/sum(df_psampmean$p) )
#実験毎のyをサンプリング
y_rep<- sample(c(1, 0), size = n[j], prob = c(theta_rep, 1-theta_rep) , replace = T)
#実験ごとのy/nを求める
yn_rep[j] <- sum(y_rep)/n[j]
}
#61番目と6番目
t_st$t_rep[i] <- n_order(yn_rep, 61) - n_order(yn_rep, 6)
#平均
t_st$t_mean[i] <- mean(yn_rep)
}
hist(yn_obs)
hist(yn_rep)
ggplot(data = t_st,aes(x = t_rep)) + geom_histogram()+geom_vline(xintercept= n_order(yn_obs, 61) - n_order(yn_obs, 6))+ggtitle("The difference between 61st and 6th of y/n")
ggplot(data = t_st,aes(x = t_mean)) + geom_histogram()+geom_vline(xintercept= mean(yn_obs))+ggtitle("The mean of y/n")
theta_rep
k <- t_st %>% filter(t_mean > mean(yn_obs)) %>% nrow()
k/200
```
結論。シミュレーションによって$y/n$の平均に近い値を出すことはできる。しかし、裾の広さは、今回のシミュレーションでは表現することができなかった。また、平均を試験量とするp値は0.6~0.7あり、サンプリングによる誤差であるとは言えない。そのため、今回の試行においては十分にデータの背景にあるモデルを特定できなかったと結論する。
改善の方向性としては、他の$\alpha, \beta$についても検証をすること、そして、$(p(\alpha, \beta)$など)モデリングそのものを検証することも考えられる。しかし、$\alpha, \beta$の複数の検証はチェリーピッキングになりかねず、また、$p(\alpha, \beta)$の他の仮定はあまり現実的とは言えない。では何がゴール?
Sign in with Wallet
Connect another wallet