\DeclareMathOperator*{\E}{\mathbb{E}}

# 5章　章末問題

## 1:

Exchangeability with known model parameters: For each of the following three examples,
answer: (i) Are observations $y_{1}$ and $y_{2}$ exchangeable? (ii) Are observations $y_{1}$ and $y_{2}$
independent? (iii) Can we act as if the two observations are independent?

### a

A box has one black ball and one white ball. We pick a ball $y_{1}$ at random, put it back,
and pick another ball $y_{2}$ at random.

#### solution

This is a problem of sampling with replacement. The observations are therefore a pair of bernoulli draws with a shared parameter $p$.

\begin{align}
y_{1} \sim Bernoulli(p) \\
y_{2} \sim Bernoulli(p)
\end{align}

The joint distribution of both observations is

\begin{cases}
(1-p)^{2} &\text{for $y_{1} = 0, y_{2} = 0$}\\
p(1-p) &\text{for $y_{1} = 1, y_{2} = 0$}\\
(1-p)p &\text{for $y_{1} = 0, y_{2} = 1$}\\
p^{2} &\text{for $y_{1} = 1, y_{2} = 1$}\\
0&\text{otherwise}\\
\end{cases}

This is equivalent to

\begin{align}
(p^{y_{1}}(1-p)^{1-y_{1}})(p^{y_{2}}(1-p)^{1-y_{2}})
\end{align}

and the observations are independent. Independence implies exchangeability so the observations are therefore exchangeable.

### a

A box has one black ball and one white ball. We pick a ball $y_{1}$ at random, we do not
put it back, then we pick ball $y_{2}$.

#### solution

The distribution of $y_{1}$ is a bernoulli draw with probability $p=k/n = 0.5$ where $k$ is the number of white balls and $n$ the total number of balls.

\begin{align}
p(y_{1}) = 0.5^{y_{1}}0.5^{1-{y_{1}}}
\end{align}

$y_{2}$ is a bernoulli draw with the probability $p = k - y_{1}/(n-1)$

\begin{align}
p(y_{2}) = (\frac{k - y_{1}}{n-1})^{y_{2}}(1 - (\frac{k - y_{1}}{n-1}))^{1 - y_{2}}
\end{align}

## 2: Exchangeability with unknown model parameters
### Question
For each of the following three exam- ples, answer: (i) Are observations $y_1$ and $y_2$ exchangeable? (ii) Are observations $y_1$ and $y_2$ independent? (iii) Can we act as if the two observations are independent?
#### (a)
A box has n black and white balls but we do not know how many of each color. We pick a ball $y_1$ at random, put it back, and pick another ball $y_2$ at random.
#### Solution
(i)yes, (ii)yes, (iii)yes
毎回ボールを戻しているので、試行毎での情報は変化しない

#### (b)
A box has n black and white balls but we do not know how many of each color. We pick a ball $y_1$ at random, we do not put it back, then we pick ball $y_2$ at random.
#### Solution
(i)no, (ii)no, (iii)no
取り出したボールを元に戻さないので、試行毎で取り出す際に知りうる情報が変化しており、exchangeableでないし、先に取り出したボールの色によって次の試行についての情報が変化するので独立した試行とも扱えない。

#### \(c\)
Same as (b) but we know that there are many balls of each color in the box.
#### Solution
(i)no, (ii)no, (iii)yes
解析的には\(c\)と同様だが、nが非常に大きい場合には独立と扱える。

## 3
#### Question
Hierarchical models and multiple comparisons:
(a) Reproduce the computations in Section 5.5 for the educational testing example. Use
the posterior simulations to estimate (i) for each school j, the probability that its
coaching program is the best of the eight; and (ii) for each pair of schools, j and k,
the probability that the coaching program in school j is better than that in school k.
(b) Repeat (a), but for the simpler model with * set to ' (that is, separate estimation
for the eight schools). In this case, the probabilities (ii) can be computed analytically.
\(c\) Discuss how the answers in (a) and (b) differ.
(d) In the model with * set to 0, the probabilities (i) and (ii) have degenerate values; what
are they?

#### Solution
##### (a)
(p.118より）学校ごとのパラメータ$\theta$の推定をするにあたっては、シミュレーションを用いる。$(\theta, \mu, \tau)$の同時事後分布は、
\begin{align}
p(\theta, \mu, \tau | y) = p(\tau | y)p(\mu, \tau | y)p(\theta| \mu, \tau, y)
\end{align}
なので、まず$\tau$の事後分布$p(\tau | y)$から$\tau$のありうる値を考え、その値を仮定した上で$\mu$の事後分布$p(\mu, \tau | y)$から$\mu$をシミュレーションし、その結果に基づいて$p(\theta| \mu, \tau, y)$から$\theta$を計算する。
###### $\tau$の推定
　まず、(5.21)より$\tau$の事後分布$p(\tau | y)$は、
\begin{align}
p(\tau| y) \propto p(\tau)\prod_{j=1}^{J}(\sigma_{j}^2 + \tau^2)^{-1/2}\exp(-\frac{(\bar{y}_{\cdot j}-\hat{\mu})^{2}}{2(\sigma_{j}^2 + \tau^2)})
\end{align}
ここで(5.20)より、
\begin{align}
\hat{\mu}=\frac{\sum_{j=1}^{J}\frac{1}{\sigma_{j}^2 + \tau^2}\bar{y}_{\cdot j}}{\sum_{j=1}^{J}\frac{1}{\sigma_{j}^2 + \tau^2}}
\end{align}

\begin{align}
V_{\mu}^{-1} = \sum_{j=1}^{J}\frac{1}{\sigma_{j}^2 + \tau^2}
\end{align}
これらを計算するRコードは以下の通り。(なお事前分布は本文通り$p(\tau)\propto 1$を仮定)

# Mu_hat
mu_hat <- function(sigma, y, tau){
  mu_vector <- NULL
  for(i in 1:length(tau)){
    tau_i     <- tau[i]
    mu        <- sum(y/((sigma)^2+(tau_i)^2))/sum(1/((sigma)^2+(tau_i)^2))
    mu_vector <- append(mu_vector, mu)
  }
  return(mu_vector)
}

# inverse of V_mu
inv_Vmu <- function(sigma, tau){
  prec_vector <- NULL
  for(i in 1:length(tau)){
    tau_i       <- tau[i]
    prec        <- sum(1/((sigma)^2 + (tau_i)^2))
    prec_vector <- append(prec_vector, prec)
  }
  return(prec_vector)
}

# Posterior of Tau
post <- function(sigma, y, tau){
  post_tau_vector <- NULL
  for (i in 1:length(tau)) {
    tau_i    <- tau[i]
    mu_hat1  <- mu_hat(sigma, y, tau_i)
    sqrt_Vmu <- inv_Vmu(sigma, tau_i)
    mu       <- sum(y/((sigma)^2+(tau_i)^2))/sum(1/((sigma)^2+(tau_i)^2))
    outer <- ((sigma)^2 + (tau_i)^2)^(-1/2)
    expon <- exp((-(y - mu_hat1)^2)/(2*(sigma^2 + tau_i^2)))
    post_tau <- sqrt_Vmu * prod(outer * expon)
    post_tau_vector <- append(post_tau_vector, post_tau)
  }
  return(post_tau_vector)
}

$\tau$の推定を行っているFigure 5.5を再現する。
（[5.21]は$\tau$について非常に複雑な形の関数となっているので、推定と言えるのかわかりませんが・・・）

# データ
y            <- c(28,  8, -3,  7, -1,  1, 18, 12) ## Y（データ）
sigma        <- c(15, 10, 16, 11,  9, 11, 10, 18) ## Sigma（データ）
tau          <- seq(0, 30, 0.01) ## Tau（変数）

# グラフ
post_tau_sim <- post(sigma, y, tau)
plot (tau, post_tau_sim, ylim=c(0,1.1*max(post_tau_sim)),
      type="l", xlab="tau", ylab="", xaxs="i",
      yaxs="i", yaxt="n", bty="n", cex=2)

本文でも言及されている通り、$\tau=0$あたりが最頻値となる。
###### $\mu$と$\theta$の推定
$\tau$、$\mu$を前提とした$\mu$の事後分布は$N(\hat{\mu}, V_{\mu})$である。ここでは$\mu$について平均を取る、すなわち$\mu=\hat{\mu}$とし、その上で$\theta$のシミュレーションを行う。
　$\theta = (\theta_1, \theta_2, ..., \theta_j , ... \theta_J)$であり、$\theta_j$は独立・交換可能である。$\theta_j$の条件づき事後分布は、$N(\hat{\theta}_{j}, V_{j})$であるが、ここで(5.17)より、
\begin{align}
\hat{\theta_j}=\frac{\frac{1}{\sigma_{j}^2}\bar{y}_{\cdot j}+\frac{1}{\tau^2}\mu}{\frac{1}{\sigma_{j}^2}+\frac{1}{\tau^2}}
\end{align}
\begin{align}
V_{j} = \frac{1}{\frac{1}{\sigma_{j}^2}+\frac{1}{\tau^2}}
\end{align}

これに基づき、$\theta_j$の推定を行っているFigure 5.6を再現する。

theta_data <- NULL

for(i in 1:length(y)){
  theta <- NULL
  y_i     <- y[i]
  sigma_i <- sigma[i]
  for (j in 2:length(tau)) {
    theta_sim <- ((y_i/(sigma_i)^2) + (mu_hat(sigma, y, tau[j])/(tau[j])^2))/((1/(sigma_i)^2) + (1/(tau[j])^2))
    theta <- append(theta, theta_sim)
  }
  theta_data <- cbind(theta_data, theta)
}
theta_data <- rbind(theta_data[1, ], theta_data) 
# tau = 0の時、theta_hatは定義できないので、仮の値を入れている。
theta_data <- as.data.frame(theta_data)

plot (tau, theta_data[,1], ylim=c(-5,1.1*max(theta_data[,1])),
      type="l", xlab="tau", ylab="Estimated Treatment Effect", xaxs="i",
      yaxs="i", bty="n", cex=2)
par(new = TRUE)
plot (tau, theta_data[,2], ylim=c(-5,1.1*max(theta_data[,1])),
      type="l", xlab="tau", ylab="", xaxs="i",
      yaxs="i", bty="n", cex=2)
par(new = TRUE)
plot (tau, theta_data[,3], ylim=c(-5,1.1*max(theta_data[,1])),
      type="l", xlab="tau", ylab="", xaxs="i",
      yaxs="i", bty="n", cex=2)
par(new = TRUE)
plot (tau, theta_data[,4], ylim=c(-5,1.1*max(theta_data[,1])),
      type="l", xlab="tau", ylab="", xaxs="i",
      yaxs="i", bty="n", cex=2)
par(new = TRUE)
plot (tau, theta_data[,5], ylim=c(-5,1.1*max(theta_data[,1])),
      type="l", xlab="tau", ylab="", xaxs="i",
      yaxs="i", bty="n", cex=2)
par(new = TRUE)
plot (tau, theta_data[,6], ylim=c(-5,1.1*max(theta_data[,1])),
      type="l", xlab="tau", ylab="", xaxs="i",
      yaxs="i", bty="n", cex=2)
par(new = TRUE)
plot (tau, theta_data[,7], ylim=c(-5,1.1*max(theta_data[,1])),
      type="l", xlab="tau", ylab="", xaxs="i",
      yaxs="i", bty="n", cex=2)
par(new = TRUE)
plot (tau, theta_data[,8], ylim=c(-5,1.1*max(theta_data[,1])),
      type="l", xlab="tau", ylab="", xaxs="i",
      yaxs="i", bty="n", cex=2)


(Figure 5.7の再現も試みましたが、うまくいきませんでした。以下没コード)

theta_sd_data <- NULL

for(i in 1:length(y)){
  theta_sd_vector <- NULL
  y_i     <- y[i]
  sigma_i <- sigma[i]
  for (j in 2:length(tau)) {
    mu_sim    <- rnorm(1000, mu_hat(sigma_i, y_i, tau[j]), (inv_Vmu(sigma_i, tau[j]))^(-1))
    theta_sim <- ((y_i/(sigma_i)^2) + (mu_sim/(tau[j])^2))/((1/(sigma_i)^2) + (1/(tau[j])^2))
    theta_sd  <- sd(theta_sim)
    theta_sd_vector <- append(theta_sd_vector, theta_sd)
  }
  theta_sd_data <- cbind(theta_sd_data, theta_sd_vector)
}
View(theta_sd_data)

###### グループ間比較
$\tau = 6$とし、グループ間比較を行う。
（教科書では$\tau$の値をどうしているか全く説明がありませんでしたが、6にすると解答に近い値が出ました。）
それぞれのグループにおける$\theta$が最も大きい確率を計算する。

set.seed(19940309)
tau_set <- 6
theta_sim_data <- NULL

# Theta_jをN(theta_hat_j, V_j)から1000ずつ発生させる。
for (i in 1:length(sigma)) {
  theta_hat <- ((y[i]/(sigma[i])^2) + (mu_hat(sigma, y, tau_set)/(tau_set)^2))/((1/(sigma[i])^2) + (1/(tau_set)^2))
  v         <- 1/((1/(sigma[i])^2) + (1/(tau_set)^2))
  theta_sim <- rnorm(1000, theta_hat, sqrt(v))
  theta_sim_data <- cbind(theta_sim_data, theta_sim)
}
colnames(theta_sim_data) <- LETTERS[1:8]
theta_sim_data <- as.data.frame(theta_sim_data)
theta_sim_data$max <- NA

# それぞれのグループが最大値になっているかどうかを識別
for (i in 1:nrow(theta_sim_data)) {
  theta_sim_data$max[i] <- max(theta_sim_data[i, ], na.rm = TRUE)
}
theta_sim_data$maxA <- ifelse(theta_sim_data$A == theta_sim_data$max, 1, 0)
theta_sim_data$maxB <- ifelse(theta_sim_data$B == theta_sim_data$max, 1, 0)
theta_sim_data$maxC <- ifelse(theta_sim_data$C == theta_sim_data$max, 1, 0)
theta_sim_data$maxD <- ifelse(theta_sim_data$D == theta_sim_data$max, 1, 0)
theta_sim_data$maxE <- ifelse(theta_sim_data$E == theta_sim_data$max, 1, 0)
theta_sim_data$maxF <- ifelse(theta_sim_data$F == theta_sim_data$max, 1, 0)
theta_sim_data$maxG <- ifelse(theta_sim_data$G == theta_sim_data$max, 1, 0)
theta_sim_data$maxH <- ifelse(theta_sim_data$H == theta_sim_data$max, 1, 0)

#　最大値になっている確率を計算
mean(theta_sim_data$maxA) # 0.274
mean(theta_sim_data$maxB) # 0.098
mean(theta_sim_data$maxC) # 0.084
mean(theta_sim_data$maxD) # 0.104
mean(theta_sim_data$maxE) # 0.048
mean(theta_sim_data$maxF) # 0.049
mean(theta_sim_data$maxG) # 0.218
mean(theta_sim_data$maxH) # 0.125

（二群間比較は時間がないので省略しました。すみません・・・）

##### (b)
$\tau\rightarrow \infty$とすると、全てのグループを別々に推定するのと同値となる。この時、事後分布は$N(y_j, \sigma_j^{2})$となる。
この時、各グループが最大になる確率は以下の通り。

theta_sim_data <- NULL
# Theta_jをN(y_j, sigma_j^2)から1000ずつ発生させる。
for (i in 1:length(y)) {
  theta_sim <- rnorm(1000, y[i], sigma[i])
  theta_sim_data <- cbind(theta_sim_data, theta_sim)
}
theta_sim_data <- as.data.frame(theta_sim_data)
colnames(theta_sim_data) <- LETTERS[1:8]

theta_sim_data$max <- NA
for (i in 1:nrow(theta_sim_data)) {
  theta_sim_data$max[i] <- max(theta_sim_data[i, ], na.rm = TRUE)
}
theta_sim_data$maxA <- ifelse(theta_sim_data$A == theta_sim_data$max, 1, 0)
theta_sim_data$maxB <- ifelse(theta_sim_data$B == theta_sim_data$max, 1, 0)
theta_sim_data$maxC <- ifelse(theta_sim_data$C == theta_sim_data$max, 1, 0)
theta_sim_data$maxD <- ifelse(theta_sim_data$D == theta_sim_data$max, 1, 0)
theta_sim_data$maxE <- ifelse(theta_sim_data$E == theta_sim_data$max, 1, 0)
theta_sim_data$maxF <- ifelse(theta_sim_data$F == theta_sim_data$max, 1, 0)
theta_sim_data$maxG <- ifelse(theta_sim_data$G == theta_sim_data$max, 1, 0)
theta_sim_data$maxH <- ifelse(theta_sim_data$H == theta_sim_data$max, 1, 0)

mean(theta_sim_data$maxA) # 0.553
mean(theta_sim_data$maxB) # 0.028
mean(theta_sim_data$maxC) # 0.032
mean(theta_sim_data$maxD) # 0.036
mean(theta_sim_data$maxE) # 0.001
mean(theta_sim_data$maxF) # 0.01
mean(theta_sim_data$maxG) # 0.167
mean(theta_sim_data$maxH) # 0.173

（ここでも二群間比較は時間がないので省略しました。すみません・・・）

##### \(c\)
$\tau \rightarrow \infty$とすると、より極端な確率が出やすい。例えばAが最大になる確率は$\tau = 6$の時0.274であるが、$\tau \rightarrow \infty$とすれば0.553にまで跳ね上がる。逆にEが最大になる確率は前者では0.048だが、後者では0.001まで下がる。そもそも$\tau$とは、$\theta_j$を正規分布から引き出された確率変数と見做した時のその正規分布の標準偏差の値である（$\theta_j \sim N(\mu, \tau^2)$、5.14参照）。したがって$\tau$を大きい値に設定すれば、$\theta_j$のばらつきが大きくなるのも当然である。

##### \(d\)
$\tau = 0$とすれば、$\theta_j$全て同じ値になる。$\theta_j \sim N(\mu, \tau^2)$に鑑みれば、この結論も当然である。

## 4
### Question
Exchangeable prior distributions: suppose it is known a priori that the $2J$ parameters $\theta_1$, . . . , $\theta_{2J}$ are clustered into two groups, with exactly half being drawn from a $N(1, 1)$ distribution, and the other half being drawn from a $N(−1, 1)$ distribution, but we have not observed which parameters come from which distribution.
(a) Are $\theta_1, . . . , \theta_{2J}$ exchangeable under this prior distribution?
(b) Show that this distribution cannot be written as a mixture of independent and iden-tically distributed components.
(c) Why can we not simply take the limit as $J \rightarrow \infty$ and get a counterexample to de Finetti's theorem?
See Exercise 8.10 for a related problem.
### Solution
#### a
交換可能性を満たしている。$\theta_1, . . . , \theta_{2J}$を区別する情報が他にないため。

#### b
5.9.5より、$\theta_1, . . . , \theta_{2J}$のうちある二つの共分散が負であることを示せば、対偶によりi.i.d.の混合分布の形式で書き表すことができないことを示せる。

今、$\theta_1, . . . , \theta_{2J}$のうち最大値を$\theta_{M}$とおく。この時、$\theta_{M} \sim N(1, 1)$である可能性が高く、それ以外の$\theta_i$は、$\theta_i \sim N(-1, 1)$である可能性が高くなる。よって、$cov(\theta_M, \theta_i)<0$であり、i.i.dではないことは示せる。

#### c
de Finetti's theorem (p. 105)とは、$J \rightarrow \infty$の時(5.2)のi.i.d式で表せることを示す命題である。すなわち、この問題は、i.i.d式の反例、すなわち$cov(\theta_i, \theta_j) < 0$にならないことを示せば良い。

$J \rightarrow \infty$の時、ある$\theta_i \sim N(1, 1)$として、別の$\theta_j \sim N(-1, 1)$が満たされる確率は変わらない。すなわち、$cov(\theta_i, \theta_j) = 0$である。よって、de Finetti's theoremの反例をとることはできない。

## 5.

式 1.9 (p.21) より、
\begin{align}
\mathrm{cov}(\theta_i, \theta_j) &= \mathrm{E}(\mathrm{cov}(\theta_i, \theta_j|\phi)) + \mathrm{cov}(\mathrm{E}(\theta_i|\phi), \mathrm{E}(\theta_j|\phi)) \\
&= 0 + \mathrm{var}(\mathrm{E}(\theta_j|\phi)) \ \ (\because \theta_j|\phi \ (j \in 1, 2, \cdots, J) \ \mathrm{is \ i.i.d.}) \\
&\geq 0.
\end{align}

## 6
### Question
Exchangeable models:
(a) In the divorce rate example of Section 5.2, set up a prior distribution for the values $y_1, . . . , y_8$ that allows for one low value (Utah) and one high value (Nevada), with independent and identical distributions for the other six values. This prior distribution should be exchangeable, because it is not known which of the eight states correspond to Utah and Nevada.
(b) Determine the posterior distribution for $y_8$ under this model given the observed values of $y_1, . . . , y_7$ given in the example. This posterior distribution should probably have two or three modes, corresponding to the possibilities that the missing state is Utah, Nevada, or one of the other six.
(c\) Now consider the entire set of eight data points, including the value for $y_8$ given at the end of the example. Are these data consistent with the prior distribution you gave in part (a) above? In particular, did your prior distribution allow for the possibility that the actual data have an outlier (Nevada) at the high end, but no outlier at the low end?

### Solution
#### (a)
事前分布は１つのlow valueと１つのhigh valueがあり、残りの州に関しては独立に同一の分布に従う。この時、exchangeabilityが満たされなければならないので、どの州がどの番号に該当するのかがわからない。なので、1つのlow及び1つのhigh valueが許されるような全ての事前分布の平均を考える。８つの州からhigh, lowをそれぞれ選ぶので、
\begin{align}
p(\theta_1, \theta_2, ...,\theta_8) = \frac{1}{56}\sum_{f}\mathrm{Beta}(\theta_{f(1)}|\alpha_l, \beta_l)\mathrm{Beta}(\theta_{f(2)}|\alpha_h, \beta_h)\prod_{i=3}^8\mathrm{Beta}(\theta_{f(i)}|\alpha_m, \beta_m)
\end{align}
ここで、$f:(1,..., 8) \rightarrow (1,..., 8)$。また、$\alpha, \beta$はハイパーパラメータ。

#### (b)
\begin{align}
p(\theta_8|\theta_1,...,\theta_7) = \frac{p(\theta_1,...,\theta_7, \theta_8)}{p(\theta_1,...,\theta_7)}
\end{align}
(a)の結果をもとにそれぞれ値を計算する。
残りはRmarkdown上で。
#### \(c\)
(事前分布のパラメータにもよるとは思いますが)分布が重なっている限り、外れ値と予想される州が外れ値でないことを事前分布は許容しているといえる。

## 7

If $y|\theta \sim Poisson(\theta)$, and $\theta \sim Gamma(\alpha, \beta)$, then the marginal (prior predictive)
distribution of $y$ is negative binomial with parameters $\alpha$ and $\beta$ (or $p = \beta)/(1 + \beta))$).
Use the formulas (2.7) and (2.8) to derive the mean and variance of the negative binomial.

- 2.7 $E(y) = E(E(y|\theta))$
- 2.8 $var(y) = E(var(y|\theta)) + var(E(y|\theta))$

#### solution

first the condiditional expectation of y given theta.

\begin{align}
E(y|\theta) = E(Poisson(\theta)) = \theta
\end{align}

We know that this is the case from the properties of the poisson distribution. This means that $E(y) = E(\theta)$. $E(\theta)$ is equal to $E(Gamma(\alpha, \beta))$

\begin{align}
E(\theta) = E(Gamma(\alpha, \beta)) = \frac{\alpha}{\beta}
\end{align}

Next, the variance. The first terms is:

\begin{align}
E(var(y|\theta)) = E(\theta) = \frac{\alpha}{\beta}
\end{align}

Because the expectation of the poisson distribution is $\theta$ and theta comes from a gamma distribution

\begin{align}
var(E(y|\theta)) = var(\theta) = \frac{\alpha}{\beta^{2}}
\end{align}

The variance term is then equal to:

\begin{align}
E(var(y|\theta)) + var(E(y|\theta)) = \frac{\alpha}{\beta} + \frac{\alpha}{\beta^{2}}
\end{align}

This can be simplified to:

\begin{align}
\frac{\alpha}{\beta} + \frac{\alpha}{\beta^{2}} = \frac{\alpha\beta + \alpha}{\beta^{2}} = \frac{\alpha}{\beta^{2}(\beta + 1)}
\end{align}

#### b solution

In the normal model with unknown location and scale $(\mu, \sigma^{2}), the noninformative
prior density, $p(\mu, \sigma^{2}) \propto 1/\sigma^{2}$, results in a normal-inverse-$\chi_{2}$ posterior distribution for
$(\mu, \sigma^{2})$. Marginally then $\sqrt n(\mu −\bar{y})/s$ has a posterior distribution that is $t_{n}−1$. Use
(2.7) and (2.8) to derive the first two moments of the latter distribution.

We are interested in knowing the
\begin{align}
$\sqrt n(\mu − \bar{y})/s$

## 8. Discrete mixture models
### Question
if $p_m(θ)$, for $m = 1,...,M$, are conjugate prior densities for the sampling model y|θ, show that the class of finite mixture prior densities given by
\begin{align}
p(\theta)=\Sigma_{m=1}^{M}\lambda_m p_m(\theta)
\end{align}
is also a conjugate class, where the $\lambda_m$’s are nonnegative weights that sum to 1. This can provide a useful extension of the natural conjugate prior family to more flexible distributional forms. As an example, use the mixture form to create a bimodal prior density for a normal mean, that is thought to be near 1, with a standard deviation of 0.5, but has a small probability of being near −1, with the same standard deviation. If the variance of each observation $y_1, . . . , y_{10}$ is known to be 1, and their observed mean is $y = −0.25$, derive your posterior distribution for the mean, making a sketch of both prior and posterior densities. Be careful: the prior and posterior mixture proportions are different.

### Solution
$p_m(\theta|y)$を$p(\theta)$に対応する事後分布であるとすると、それぞれの$m$に対して
\begin{align}
p_m(\theta|y)=\frac{p_m(\theta)p(y|\theta)}{p_m(y)}
\end{align}
このとき、$p_m(y)$は事前予測分布である。
$p(\theta)=\Sigma_{m=1}^{M}\lambda_m p_m(\theta)$より、$\theta$の事後分布は、
\begin{align}
p(\theta|y) &\propto \Sigma_m \lambda_m p_m(\theta)p(y|\theta)\\
&= \Sigma_m \lambda_m p_m(y) p_m(\theta|y)
\end{align}

theta <- seq(-3,3,.01)
prior <- c (0.9, 0.1)
dens <- prior[1]*dnorm(theta,1,0.5) +
  prior[2]*dnorm(theta,-1,0.5)
plot (theta, dens, ylim=c(0,1.1*max(dens)),
      type="l", xlab="theta", ylab="", xaxs="i",
      yaxs="i", yaxt="n", bty="n", cex=2) 
mtext ("prior density", cex=2, 3)

marg <- dnorm(-.25,c(1,-1),sqrt(c(0.5,0.5)^2+1/10)) 
posterior <- prior*marg/sum(prior*marg)
dens <- posterior[1]*dnorm(theta,1.5/14,sqrt(1/14)) + 
  posterior[2]*dnorm(theta,-6.5/14,sqrt(1/14))
plot (theta, dens, ylim=c(0,1.1*max(dens)), type="l", 
      xlab="theta", ylab="", xaxs="i", yaxs="i", yaxt="n", bty="n", cex=2)
mtext ("posterior density", cex=2, 3)

![Uploading file..._educmfrzl]()

![Uploading file..._v604rg6gz]()

## 9
### Question
Noninformative hyperprior distributions: consider the hierarchical binomial model in Section 5.3. Improper posterior distributions are, in fact, a general problem with hierarchical models when a uniform prior distribution is specified for the logarithm of the population standard deviation of the exchangeable parameters. In the case of the beta population distribution, the prior variance is approximately $(\alpha+\beta)^{−1}$ (see Appendix A), and so a uniform distribution on $\log(\alpha+\beta)$ is approximately uniform on the log standard deviation. The resulting unnormalized posterior density (5.8) has an infinite integral in the limit as the population standard deviation approaches 0. We encountered the problem again in Section 5.4 for the hierarchical normal model.

#### (a)
Show that, with a uniform prior density on $(\log(\frac{\alpha}{\beta})), \log(\alpha+\beta))$, the unnormalized posterior density has an infinite integral.

#### (b)
A simple way to avoid the impropriety is to assign a uniform prior distribution to the standard deviation parameter itself, rather than its logarithm. For the beta population distribution we are considering here, this is achieved approximately by assigning a uniform prior distribution to $(\alpha+\beta)^{−1/2}$. Show that combining this with an independent uniform prior distribution on $\frac{\alpha}{\alpha+\beta}$ yields the prior density (5.10).
#### \(c\)
Show that the resulting posterior density (5.8) is proper as long as $0 < y_j < n_j$ for at least one experiment $j$.
### Solution
#### (a)
Section 5.3において$\alpha, \beta$とは、グループごとのデータ$y_j$を生み出す二項分布$Bin(n_j, \theta_j)$のパラメータ$\theta_j$が従う$Beta(\alpha, \beta)$のハイパーパラメータである。
ここで$\alpha, \beta$の事後分布は(5.8)より、
\begin{align}
p(\alpha, \beta|y) &\propto p(\alpha, \beta)\prod_{j=1}^{J}\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\frac{\Gamma(\alpha+y_j)\Gamma(\beta+n_j-y_j)}{\Gamma(\alpha+\beta+n_j)} \\
(&= p(\alpha, \beta)p(y |\alpha, \beta))
\end{align}
と表される。尤度部分$p(y |\alpha, \beta)$を見てみると、
\begin{align}
p(y |\alpha, \beta) &= \prod_{j=1}^{J}\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\frac{\Gamma(\alpha+y_j)\Gamma(\beta+n_j-y_j)}{\Gamma(\alpha+\beta+n_j)} \\
&\propto \prod_{j=1}^{J}\frac{[\alpha\cdots(\alpha+y_j-1)][\beta\cdots(\beta+n_j-y_j-1)]}{(\alpha+\beta)\cdots(\alpha+\beta+n_j-1)} \\
&\approx \prod_{j=1}^{J}\frac{\alpha^{y_j}\beta^{n_j-y_j}}{(\alpha + \beta)^{n_j}}\\
&= \prod_{j=1}^{J} \Big(\frac{\alpha}{\alpha + \beta}\Big)^{y_j} \Big(\frac{\beta}{\alpha + \beta}\Big)^{n_j-y_i}
\end{align}
ここで$y_j, n_j, \alpha/\beta$を固定すると、上記は定数となる。
したがって示すべきは、$p(\alpha, \beta)$がinfinite integralを持つことである。この問題でこれは$(\log(\frac{\alpha}{\beta})), \log(\alpha+\beta))$における一様分布なので、$\alpha + \beta \rightarrow \infty$とすればinfinite integralを持つ。（ここよくわからない？）
したがって題意は示された。

#### (b)
$(\alpha, \beta)$から$(\frac{\alpha}{\alpha+\beta}, (\alpha + \beta)^{-1/2})$の変数変換を考えると、この時のヤコビアンは、
\begin{align}
\begin{vmatrix}
\frac{\partial}{\partial\alpha}\big(\frac{\alpha}{\alpha+\beta}\big) & \frac{\partial}{\partial\beta}\big(\frac{\alpha}{\alpha+\beta}\big) \\
\frac{\partial}{\partial\alpha}\big(\alpha + \beta\big)^{-1/2} & \frac{\partial}{\partial\beta}\big(\alpha + \beta\big)^{-1/2}
\end{vmatrix} \\=
\begin{vmatrix}
\frac{\beta}{\alpha+\beta} & \frac{-\alpha}{\alpha+\beta} \\
-\frac{1}{2}(\alpha+\beta)^{-3/2} & -\frac{1}{2}(\alpha+\beta)^{-3/2}
\end{vmatrix}
\end{align}
\begin{align}
&= \Big(\frac{\beta}{\alpha+\beta}\Big)\times\Big(-\frac{1}{2}(\alpha+\beta)^{-3/2} \Big)-\Big(\frac{-\alpha}{\alpha+\beta}\Big)\times\Big(-\frac{1}{2}(\alpha+\beta)^{-3/2}\Big) \\
&= -\frac{1}{2}(\alpha+\beta)^{-5/2}
\end{align}
これで(5.9)は示された。
（ここから(5.10)にどのように行けばいいかがわかりませんでした・・・）

##### 齋藤答案
そもそも、上記答案の趣旨は$(\alpha/(\alpha + \beta), (\alpha + \beta)^{-1/2})$がそれぞれ一様分布に従うとき、$(\frac{\alpha}{\alpha+\beta}, (\alpha + \beta)^{-1/2}) \rightarrow (\alpha, \beta)$の変数変換を考えると
\begin{align}
&\int p(\alpha/(\alpha + \beta), (\alpha + \beta)^{-1/2}) d \alpha/(\alpha + \beta) d(\alpha + \beta)^{-1/2} \propto \int |J| d\alpha d\beta \propto 1 &(\because p(\alpha/(\alpha + \beta), (\alpha + \beta)^{-1/2}) \propto 1)
\end{align}
が成り立つ。よって、$p(\alpha, \beta)\propto |J|$が成り立った。そのため、ヤコビアンを求めればよかった。

同様に、(5.10)式も、$(\frac{\alpha}{\alpha+\beta}, (\alpha + \beta)^{-1/2}) \rightarrow (\log \frac{\alpha}{\beta}, \log (\alpha + \beta))$の変数変換を考え、ヤコビアンを求めれば同時分布を求められる。今、$\omega = \log \frac{\alpha}{\beta}, \gamma = \log (\alpha + \beta)$とおく。すると、$\alpha = \frac{e^{\omega + \gamma}}{e^\omega + 1}, \beta = \frac{e^\gamma}{e^\omega + 1}$である。よって、
\begin{align}
|J|&= \begin{vmatrix}
\frac{d}{d\omega}\frac{e^{\omega + \gamma}}{e^\gamma(e^\omega + 1)} &
\frac{d}{d\gamma}\frac{e^{\omega + \gamma}}{e^\gamma(e^\omega + 1}\\
\frac{d}{d\omega}e^\gamma &
\frac{d}{d\gamma}e^\gamma
\end{vmatrix}\\
&= \begin{vmatrix}
\frac{e^{\omega + 2\gamma}}{(e^\gamma(e^\omega+1 ))^2} &
0\\
0 &
\frac{-1}{2}(\frac{e^\omega +1}{e^\gamma + e^{\omega+ \gamma}})^{1/2}
\end{vmatrix}\\
&\propto \frac{e^{\omega+2\gamma}}{e^{5\gamma/2}(e^\omega + 1)^{2}}
\end{align}
そして、$\omega = \log \frac{\alpha}{\beta}, \gamma = \log (\alpha + \beta)$なので、
\begin{align}
P(\log \frac{\alpha}{\beta}, \log (\alpha + \beta)) \propto \alpha \beta(\alpha + \beta)^{-5/2}
\end{align}

#### \(c\)
(5.8)について、4つの極限を考える。$y_j>0$が成立している実験が$J_0$個、$y_j<n_j$が成立している実験が$J_1$個、$0<y_j<n_j$が成立している実験が$J_{01}$個存在するとする（？）。なお、(5.8)は以下の通り。
\begin{align}
p(\alpha, \beta|y) &\propto p(\alpha, \beta)\prod_{j=1}^{J}\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\frac{\Gamma(\alpha+y_j)\Gamma(\beta+n_j-y_j)}{\Gamma(\alpha+\beta+n_j)} \\
&\approx p(\alpha, \beta)\prod_{j=1}^{J} \Big(\frac{\alpha}{\alpha + \beta}\Big)^{y_j} \Big(\frac{\beta}{\alpha + \beta}\Big)^{n_j-y_i}
\end{align}

##### ①$\alpha+\beta$を固定し、$\alpha\rightarrow0$
この時、尤度$p(\alpha, \beta|y)$は、$\alpha^{J_0}$の部分以外定数となるので（？）、$\alpha\rightarrow0$ならば、$J_0>0$の時尤度は0となり、$J_0=0$の時定数となる。事前分布$p(\alpha, \beta)$は$\alpha\rightarrow0$ならば定数となる。したがって$J_0>0$あるいは$J_0=0$において、事後分布はfinite integral を持つ。

##### ②$\alpha+\beta$を固定し、$\beta\rightarrow0$
①と同様。

##### ③$\alpha/\beta$を固定し、$\alpha+\beta \rightarrow 0$
##### ④$\alpha/\beta$を固定し、$\alpha+\beta \rightarrow \infty$
（ここから全然わかりませんでした。）

## 10
### Question
Checking the integrability of the posterior distribution: consider the hierarchical normal model in Section 5.4.

(a) If the hyperprior distribution is $p(\mu, \tau) \propto \tau^{ −1}$ (that is, $p(\mu, \log \tau) \propto 1)$, show that the posterior density is improper.

(b) If the hyperprior distribution is $p(\mu, \tau) \propto 1$, show that the posterior density is proper if $J > 2$.

(c) How would you analyze SAT coaching data if $J = 2$ (that is, data from only two schools)?

### Solution
#### a
この問題では、$p(\theta, \mu, \tau|y)$が不適切、すなわち積分すると有限にならないことを示せば良い。今、
\begin{align}
p(\theta, \mu, \tau|y) &= \frac{p(y)p(\tau|y)p(\mu|\tau, y)p(\theta|\mu, \tau, y)}{p(y)} &(\text{Bayes' rule)}\\
&= p(\tau|y)p(\mu|\tau, y)p(\theta|\mu, \tau, y)
\end{align}
である。今、
\begin{align}
\mu|\tau, y \sim N(\hat{\mu}, V_j)
\end{align}
ただし、
\begin{align}
\hat{\mu} = \frac{\sum_j^J \frac{1}{\sigma_j^2+\tau^2}\bar{y}_{.j}}{\sum_j^J \frac{1}{\sigma_j^2+\tau^2}}
&\text{ and }&V_{\mu}^{-1} = \sum_j^J \frac{1}{\sigma_j^2+\tau^2} \tag{5. 20}
\end{align}
である。そして、
\begin{align}
\theta_j|\mu, \tau, y \sim N(\hat{\theta_j}, V_j)
\end{align}
ただし、
\begin{align}
\hat{\theta_j} = \frac{\frac{1}{\sigma_j^2}\bar{y}_{.j}+\frac{1}{\tau^2}\mu}{\frac{1}{\sigma_j^2}+\frac{1}{\tau^2}} \tag{5.17}
\end{align}
すなわち、$p(\mu|\tau, y), p(\theta|\mu, \tau, y)$はいずれも適切な確率分布である。よって、$p(\tau|y)$が適切/不適切な確率分布であれば、$p(\theta, \mu, \tau|y)$も適切/不適切である。

今、
\begin{align}
p(\tau|y) \propto p(\tau)V_\mu^{1/2} \prod_{j = 1} (\sigma_j^2+\tau^2)^{-1/2} \exp(-\frac{(\bar{y}_{.j}-\hat{\mu})^2}{2(\sigma_j^2 + \tau^2)}) \tag{5.21}
\end{align}
である。

そして、$p(\tau) \propto p(\mu, \tau) \propto \tau^{-1}$である。(5.21)式の$p(\tau)$を除く部分について、0への極限を取ると、$\tau$について定数になり、また、$\tau^{-1}$は0への極限において無限に発散するので、$\lim_{\tau \rightarrow 0} p(\tau|y)$は発散する。よって、$p(\tau|y)$は不適切な確率分布であり、$p(\theta, \mu, \tau|y)$もまた不適切である。

#### b
aと同様に、$p(\tau|y)$が適切な確率分布かどうかを検証する。

今、$p(\tau) \propto p(\mu, \tau) \propto 1$より、(5.21)式は0近傍について発散しない。次に、$\tau \rightarrow \infty$を検証する。まず、$\lim_{\tau \rightarrow \infty} \exp(-\frac{(\bar{y}_{.j}-\hat{\mu})^2}{2(\sigma_j^2 + \tau^2)}) = 1$である。よって、$V_\mu^{1/2} \prod_{j = 1} (\sigma_j^2+\tau^2)^{-1/2}$について考えれば良い。

\begin{align}
\lim_{\tau \rightarrow \infty} p(\tau|y) &\propto p(\tau)(\sum_j \frac{1}{\sigma_j^2 + \tau^2} \prod_{j = 1} (\sigma_j^2+\tau^2))^{-1/2} \\
&= p(\tau) ( \frac{1}{\sigma_1^2 + \tau^2} (\sigma_1^2+\tau^2) (\sigma_2^2+\tau^2)\cdots+\frac{1}{\sigma_2^2 + \tau^2} (\sigma_1^2+\tau^2) (\sigma_2^2+\tau^2)\cdots+\cdots)^{-1/2}\\
&= p(\tau) (\sum_j\prod_{k\neq j}(\sigma_k^2+\tau^2))^{-1/2}\\
&\leq p(\tau)(\sum_j\prod_{k\neq j}(\tau^2))^{-1/2} &(\sigma^2 \geq 0)\\
&= \frac{p(\tau)}{J^{1/2}\tau^{J-1}} \\
&\propto \frac{1}{J^{1/2}\tau^{J-1}} &(p(\tau) \propto 1)
\end{align}

今、$J=1$の時は、$\lim_{\tau \rightarrow \infty} p(\tau|y) \propto 1$であり、0にならず、積分は発散する。$J > 1$の時は、$\lim_{\tau \rightarrow \infty} \frac{1}{J^{1/2}\tau^{J-1}} =0$であり、積分可能である。

#### c
（以下、解答DeepLコピペ)
ここではいくつかの合理的な選択肢があります。1つのアプローチは、階層モデルをやめて、2つの学校に独立した無情報事前分布を当てはめることです（練習問題3.3のように、ただし分散は既知）。
もう1つのアプローチは、τの一様事前分布を続け、結果として生じる不適切な事後分布の難しさを解析的に解決しようとすることです。この場合、事後分布は密度がすべてτ = ∞の近くにあるので、解析ではShrinkageが起こりません（(5.17)と(5.20)参照）。この分析は、実際には、学校を独立して分析する最初のアプローチと同じです。
無情報事前分布に基づく分析が十分に正確でない場合、他の学校での以前のコーチング実験の分析などの外部の知識に基づいて、τ（あるいは（µ, τ））に有情報事前分布を割り当てることが価値あることかもしれません。しかし、ここで注意しなければならないのは、2つの学校のデータしかない場合、推論は、手元のデータからは確認しにくい事前の仮定の影響を受けやすくになるということです。

## 11.

### a.

\begin{align}
p(\theta, \mu, \tau|y) \propto p(y|\theta)p(\theta|\mu, \tau)p(\mu, \tau).
\end{align}

ここで、
\begin{align}
p(y|\theta) \propto \prod_{j = 1}^J \theta_j^{y_j} (1 - \theta_j)^{n_j - y_j}.
\end{align}

また、$p(\theta|\mu, \tau)$は、変数変換の公式を用いることにより、
\begin{align}
\frac{d\mathrm{logit}(\theta_j)}{d\theta} &= \left(\log\left(\frac{\theta_j}{1 - \theta_j}\right)\right)' \\
&= \frac{1 - \theta_j}{\theta_j} \times \frac{(1 - \theta_j)\theta_j' - \theta_j(1 - \theta_j)'}{(1 - \theta_j)^2} \\
&= \theta_j^{-1} (1 - \theta_j)^{-1}
\end{align}
より、
\begin{align}
p(\theta|\mu, \tau) \propto \prod_{j = 1}^J \theta_j^{-1} (1 - \theta_j)^{-1}\tau^{-1} \exp\left(-\frac{1}{2}\left(\frac{\mathrm{logit}(\theta_j) - \mu}{\tau}\right)^2\right)
\end{align}

## 12
### Question
Conditional posterior means and variances: derive analytic expressions for $E(\theta_j |\tau, y)$ and var$(\theta_j |\tau, y)$ in the hierarchical normal model (and used in Figures 5.6 and 5.7). (Hint: use (2.7) and (2.8), averaging over $\mu$.)

### Solution
(2.7)はLaw of Iterated Expectation, (2.8)はLaw of Total Variance。これらを用いると
\begin{align}
E[\theta_j|\tau,y] &= E[E[\theta_j|\mu,\tau,y]|\tau,y]\\
&= E[\frac{y_j/\sigma_j^2 + \mu/\tau^2}{1/\sigma_j^2+1/\tau^2}|\tau, y]\\
&= \frac{y_j/\sigma_j^2 + \hat\mu/\tau^2}{1/\sigma_j^2+1/\tau^2}\\
\end{align}
なお、１行目は(5.11)、２行目は(5.17)参照。ここでは、分散既知、つまり$\sigma_j$はわかっています(p.113参照)。

\begin{align}
var(\theta_j|\tau,y) &= E[var(\theta_j|\mu,\tau,y)|\tau, y] + var[E(\theta_j|\mu,\tau,y)|\tau, y)]\\
&= E\biggl[\frac{1}{1/\sigma^2 + 1/\tau^2}\biggr] + V\biggl[\frac{y_j/\sigma_j^2 + \mu/\tau^2}{1/\sigma_j^2+1/\tau^2} \biggr]\\
&=E\biggl[\frac{1}{1/\sigma^2 + 1/\tau^2}\biggr] + V\biggl[\frac{\mu/\tau^2}{1/\sigma_j^2+1/\tau^2} \biggr]\\
&=E\biggl[\frac{1}{1/\sigma^2 + 1/\tau^2}\biggr] + (\frac{1/\tau^2}{1/\sigma^2 + 1/\tau^2})^2V(\mu)\\
&= \frac{1}{1/\sigma^2 + 1/\tau^2} +(\frac{1/\tau^2}{1/\sigma^2 + 1/\tau^2})^2 V_\mu
\end{align}
第一項も第二項も(5.17)を代入しただけ。

## 14. Hierarchical Poisson model

## 14. Hierarchical Poisson model
### Question
consider the dataset in the previous problem, but suppose only the total amount of traffic at each location is observed.

#### (a)
Set up a model in which the total number of vehicles observed at each location j follows a Poisson distribution with parameter $θ_j$, the ‘true’ rate of traffic per hour at that location. Assign a gamma population distribution for the parameters $θ_j$ and a noninformative hyperprior distribution. Write down the joint posterior distribution.

#### Solution

\begin{aligned}
p(\theta,\alpha,\beta|y) &\propto p(\alpha, \beta) \times p(\theta|\alpha, \beta) \times p(y|\theta,\alpha, \beta) \nonumber \\
p(\theta,\alpha,\beta|y) &\propto 1\times \prod_{j=1}^{10}Gamma(\theta_{j} | \alpha, \beta) \times \prod_{j=1}^{10}Poisson(y_{j}|\theta_{j}) \nonumber \\
&= \prod_{j=1}^{10}\frac{\beta^{\alpha}}{\Gamma(\alpha)}\theta_{j}^{\alpha-1}exp(-\beta \theta) \times \frac{\theta_{j}^{y_{i}}exp(-\theta_{j})}{!y_{j}} \nonumber \\
&\propto \frac{\beta^{n\alpha}}{\Gamma(\alpha)^{n}}exp(-\sum \theta_{j}( 1 + \beta )) \prod_{j=1}^{10} \theta_{j}^{\alpha + y_{j}-1} \label{bic.joint.post2}
\end{aligned}

#### (b)
Compute the marginal posterior density of the hyperparameters and plot its contours. Simulate random draws from the posterior distribution of the hyperparameters and make a scatterplot of the simulation draws.
#### Solution
グリッドで描画するために、$p(\alpha, \beta|y)$がわかれば良いので、$p(\theta, \alpha, \beta|y)=p(\alpha, \beta|y)p(\theta|\alpha, \beta, y)$から、$p(\theta|\alpha, \beta, y)$と$p(\theta, \alpha, \beta|y)$の比がわかれば良い。

Then compute the marginal posterior of $\theta$, conditional on $\alpha, \beta$. For the gamma-poisson case we have that given the hyper-parameters, each $\theta_{j}$ has a posterior distribution $Gamma(\alpha + n_{j}, \beta +1)$. Assuming exchangeability:

\begin{align}
p(\theta|\alpha,\beta,y) &\propto \prod_{j=1}^{10}Gamma(\theta_{j} | \alpha +y_{j}, \beta+1) \nonumber \\
&\propto \prod_{j=1}^{10} \theta_{j}^{\alpha + y_{j} -1}exp(-(\beta+1) \theta_{j})
\label{bic.cond.post.theta2}
\end{align}

(a)の結果と合わせると、

$$
\begin{aligned}
p(\alpha,\beta|y) &\propto \frac{\beta^{n\alpha}}{\Gamma(\alpha)^{n}} \prod_{i=1}^{n}\frac{\Gamma(\alpha+y_{i})}{(\beta + 1)^{\alpha+y_{i}}}
\label{bic.marg.post.phi}
\end{aligned}
$$

ここで、$\alpha$と$\beta$のおおよその範囲を確定させるためにデータの平均と分散から当たりをつける

$$
\begin{aligned}
\hat{\mu} &= `r mean(n)` = \frac{\hat{\alpha_{0}}}{\hat{\beta_{0}}}\\
\hat{\sigma^2} &= `r sd(n)^2` = \frac{\hat{\alpha_{0}}}{\hat{\beta_{0}}^2}
\end{aligned}
$$
Solving for $(\hat{\alpha_{0}},\hat{\beta_{0}})$:

#### \(c\)
Is the posterior density integrable? Answer analytically by examining the joint posterior density at the limits or empirically by examining the plots of the marginal posterior density above.

#### (d)
If the posterior density is not integrable, alter it and repeat the previous two steps.

#### Solution
基本線は事前分布のセッティングを変更する？

#### (e)
Draw samples from the joint posterior distribution of the parameters and hyperparameters, by analogy to the method used in the hierarchical binomial model.

#### Solution

## 15
### Question
Meta-analysis: perform the computations for the meta-analysis data of Table 5.4.
(a) Plot the posterior density of $\tau$ over an appropriate range that includes essentially all of the posterior density, analogous to Figure 5.5.
(b) Produce graphs analogous to Figures 5.6 and 5.7 to display how the posterior means and standard deviations of the $\theta_j$’s depend on $\tau$.
\(c\) Produce a scatterplot of the crude effect estimates vs. the posterior median effect
estimates of the 22 studies. Verify that the studies with smallest sample sizes are partially pooled the most toward the mean.
(d) Draw simulations from the posterior distribution of a new treatment effect, $\widetilde{\theta_j}$ . Plot
a histogram of the simulations.
(e) Given the simulations just obtained, draw simulated outcomes from replications of a hypothetical new experiment with 100 persons in each of the treated and control groups. Plot a histogram of the simulations of the crude estimated treatment effect (5.23) in the new experiment.

### Solution
Table 5.4からデータは以下のとおり。

y            <- c( 0.028, -0.741, -0.541, -0.246,  0.069,
                  -0.584, -0.512, -0.079, -0.424, -0.335,
                  -0.213, -0.039, -0.593,  0.282, -0.321,
                  -0.135,  0.141,  0.322,  0.444, -0.218,
                  -0.591, -0.608) ## Log-odds
sigma        <- c(0.850, 0.483, 0.565, 0.138, 0.281,
                  0.676, 0.139, 0.204, 0.274, 0.117,
                  0.195, 0.229, 0.425, 0.205, 0.298,
                  0.261, 0.364, 0.553, 0.717, 0.260,
                  0.257, 0.272) ## sd
tau          <- seq(0, 0.6, 0.001) ## Tau

問3と同様、関数を以下のように定義。

inv_Vmu <- function(sigma, tau){
  prec_vector <- NULL
  for(i in 1:length(tau)){
    tau_i       <- tau[i]
    prec        <- sum(1/((sigma)^2 + (tau_i)^2))
    prec_vector <- append(prec_vector, prec)
  }
  return(prec_vector)
}

mu_hat <- function(sigma, y, tau){
  mu_vector <- NULL
  for(i in 1:length(tau)){
    tau_i     <- tau[i]
    mu        <- sum(y/((sigma)^2+(tau_i)^2))/sum(1/((sigma)^2+(tau_i)^2))
    mu_vector <- append(mu_vector, mu)
  }
  return(mu_vector)
}

post <- function(sigma, y, tau){
  post_tau_vector <- NULL
  for (i in 1:length(tau)) {
    tau_i    <- tau[i]
    mu_hat1  <- mu_hat(sigma, y, tau_i)
    sqrt_Vmu <- sqrt((inv_Vmu(sigma, tau_i))^(-1))
    mu       <- sum(y/((sigma)^2+(tau_i)^2))/sum(1/((sigma)^2+(tau_i)^2))
    outer <- ((sigma)^2 + (tau_i)^2)^(-1/2)
    expon <- exp((-(y - mu_hat1)^2)/(2*(sigma^2 + tau_i^2)))
    post_tau <- sqrt_Vmu * prod(outer * expon)
    post_tau_vector <- append(post_tau_vector, post_tau)
  }
  return(post_tau_vector)
}

#### a
以下のコードでFig 5.6と同じものを再現。

post_tau_sim <- post(sigma, y, tau)
plot (tau, post_tau_sim, ylim=c(0,1.1*max(post_tau_sim)),
      type="l", xlab="tau", ylab="", xaxs="i",
      yaxs="i", yaxt="n", bty="n", cex=2)

1.3あたりが最頻値。

#### b

theta_data <- NULL

for(i in 1:length(y)){
  theta <- NULL
  y_i     <- y[i]
  sigma_i <- sigma[i]
  for (j in 2:length(tau)) {
    theta_sim <- ((y_i/(sigma_i)^2) + (mu_hat(sigma, y, tau[j])/(tau[j])^2))/((1/(sigma_i)^2) + (1/(tau[j])^2))
    theta <- append(theta, theta_sim)
  }
  theta_data <- cbind(theta_data, theta)
}
theta_data <- rbind(theta_data[1, ], theta_data) # If tau = 0, then the value of theta is not defined.
theta_data <- as.data.frame(theta_data)

frame()
plot (tau, theta_data[,1], ylim=c(-0.6,0.4),
      type="l", xlab="tau", ylab="Estimated Treatment Effect", xaxs="i",
      yaxs="i", bty="n", cex=2)
par(new = TRUE)
for (i in 2:length(y)) {
  plot (tau, theta_data[,i], ylim=c(-0.6,0.4),
        type="l", xlab="", ylab="", xaxs="i",
        yaxs="i", xaxt = "n", yaxt = "n", bty="n", cex=2)
  par(new = TRUE)
}

$\tau$が小さくなるほどグループ間の差異がなくなり、大きくなるほど差異が目立つことがわかる。
（Fig 5.7の再現はできませんでした。）

#### c

tau_set <- 0.13 # medianに設定
theta_sim_data <- NULL

for (i in 1:length(sigma)) {
  theta_hat <- ((y[i]/(sigma[i])^2) + (mu_hat(sigma, y, tau_set)/(tau_set)^2))/((1/(sigma[i])^2) + (1/(tau_set)^2))
  v         <- 1/((1/(sigma[i])^2) + (1/(tau_set)^2))
  theta_sim <- rnorm(1000, theta_hat, sqrt(v))
  theta_sim_data <- cbind(theta_sim_data, theta_sim)
}

colnames(theta_sim_data) <- paste(rep("study", length(y)), 1:length(y), sep = "")
sim_y <- apply(theta_sim_data, 2, mean)
frame()
sim_y_data <- data.frame(
  study = 1:22,
  y     = y,
  sim_y = sim_y
)
plot(sim_y_data$y, sim_y_data$sim_y,
     type = "n", xlim = c(min(y), max(y)), ylim = c(min(y), max(y)),
     xlab = "Crude Estemated Effects (Log-odds)",
     ylab = "Posterior Estimated Effect (Theta)")
text(sim_y_data$y, sim_y_data$sim_y, sim_y_data$study)

crude estimateでは大きかったばらつきが、posterior estimateでは小さくなっている（=shrinkしている）。
（Verify that the studies with smallest sample sizes are partially pooled the most toward the mean.←意味がわかりませんでした。）

#### d
ここではそれぞれのグループにおける平均パラメータ$theta_j$が共通の分布からランダムに取り出されていると考える。
この時、5.14より$\theta_j\sim N(\mu, \tau^2)$で、先ほどの分析より$\tau=1.3$、$\mu= \hat{\mu}=\frac{\sum_{j=1}^{J}\frac{1}{\sigma_{j}^2 + \tau^2}\bar{y}_{\cdot j}}{\sum_{j=1}^{J}\frac{1}{\sigma_{j}^2 + \tau^2}}$とし、$\theta_j\sim N(\hat{\mu}, 1.3^2)$から$\widetilde{\theta_j}$を1000個取り出す。（これで正しいかはわかりません・・・）

theta_pred <- rnorm(1000, mean = mu_hat(sigma, y, tau_set), sd = (inv_Vmu(sigma, tau_set))^(-1))
hist(theta_pred, freq = TRUE, 
     main = "",
     xlab = "Predicted Theta")


#### e
全然わかりませんでした・・・

## 16
### Question
Equivalent data: Suppose we wish to apply the inferences from the meta-analysis example in Section 5.6 to data on a new study with equal numbers of people in the control and treatment groups. How large would the study have to be so that the prior and data were weighted equally in the posterior inference for that study?

### Solution
p. 126以下のSAT指導の例に基づいて考える。以下は、p. 113以下に基づく。

$\sigma_j^2 = \sigma^2/n_j$として、データについて以下の分布を仮定する。
\begin{align}
y_j| \theta_j, \sigma_j \sim N(\theta_l. \sigma_j^2)
\end{align}
そして、$\theta_j$のハイパーパラメータを$\mu$、$\tau$として、分布を以下のように仮定する。
\begin{align}
\theta_j| \mu, \tau \sim N(\mu. \tau^2)
\end{align}
この時、$\bar{y}_{.j} = \frac{1}{n_j}\sum_{i=1} y_{ij}$として、事後分布を求めると、
\begin{align}
&p(\mathbf{\theta}|\mu, \tau, \mathbf{y}) \propto \prod N(\theta_j|\mu, \tau^2)\prod N(\bar{y}_{.j}|\theta_j, \sigma_j^2) \tag{5.16} \\
&\Rightarrow \theta_j|\mu, \tau, \bar{y}_{.j} \sim N(\hat{\theta}_j, V_j)
\end{align}
ただし、
\begin{align}
\hat{\theta}_j = \frac{\frac{1}{\sigma_j^2}\bar{y}_{.j}+\frac{1}{\tau^2}\mu}{\frac{1}{\sigma_j^2}+\frac{1}{\tau^2}} && V_j = \frac{1}{\frac{1}{\sigma_j^2}+\frac{1}{\tau^2}}\tag{5.17}
\end{align}

これは、p. 116の通り、事前分布$\theta_j$の平均$\mu$とj番目のグループの標本平均$\bar{y}_{.j}$とのprecision weighted平均であると言える。ここで、題意の通り、事前分布とデータとが同じ重み付けをされる条件とは、$1/\sigma_j^2 = 1/\tau^2$である。よって、$\sigma_j^2 = \sigma^2/n_j$より、
\begin{align}
n_j = \frac{\sigma^2}{\tau^2}
\end{align}