# 4章　章末問題

## 1
#### Question
Normal approximation: suppose that $y_1, ... , y_5$ are independent samples from a Cauchy distribution with unknown center $\theta$ and known scale 1: $p(y_i|\theta) \propto \frac{1}{(1 + (y_i − \theta)^2)}$. Assume that the prior distribution for $\theta$ is uniform on [0, 1]. Given the observations $(y_1,... , y_5) = (−2,−1, 0, 1.5, 2.5)$:

(a) Determine the derivative and the second derivative of the log posterior density.

(b) Find the posterior mode of $\theta$ by iteratively solving the equation determined by setting the derivative of the log-likelihood to zero.

(c\) Construct the normal approximation based on the second derivative of the log posterior density at the mode. Plot the approximate normal density and compare to the exact density as computed using the approach described in Exercise 2.11.

#### Solution
(a)
一階微分は
\begin{align}
\frac{d}{d\theta}\mathrm{log} \ p(y|\theta) &=
\frac{d}{d\theta}\mathrm{log}\frac{1}{(1 + (y_i − \theta)^2)}\\
&= \frac{2(y_i-\theta)}{1+(y_i - \theta)^2}
\end{align}

二階微分は
\begin{align}
\frac{d^2}{d\theta^2}\mathrm{log} \ p(y|\theta) &= \frac{2(\theta^2-2\theta y_i + y_i^2-1)}{(\theta^2-2\theta y_i + y_i^2+1)^2}
\end{align}

(b)
今正規分布への近似をしたいので、最頻値において一回微分が０と仮定する(p.83, ラプラス近似)。この時、分布の一回微分が0であるので、対数尤度の一回微分も０である。すなわち、
\begin{align}
\frac{d}{d\theta}\mathrm{log} \ p(y|\theta) = \sum_{i=1}^5 \frac{2(y_i-\theta)}{1+(y_i - \theta)^2} = 0
\end{align}

y <- c(-2, -1, 0, 1.5, 2.5)　#observed data
log_posterior_y <- function(y) {
  function(theta) {
    return(sum(log(1/(1 + (y - theta)^2))))
  }
} #yについての関数にする

log_posterior <- log_posterior_y(y)

mode_numerical <- optimise(log_posterior, c(0, 1), 
maximum = TRUE)$maximum

0.00006610696とほとんど0に近い値が出てくる。

(c\)
正規近似を行う。まず2階微分は
$$
\frac{2(\theta^2-2\theta y_i + y_i^2-1)}{(\theta^2-2\theta y_i + y_i^2+1)^2}
$$
であった。これに先ほどの$\hat\theta = 0$を代入すると、
$$
\sum_{i=1}^5 \frac{2(y_i^2-1)}{(y_i^2+1)^2}
$$
これが正規分布の精度パラメータ(分散の逆数)に対応する。

## 2
#### Question
Normal approximation: derive the analytic form of the information matrix and the normal approximation variance for the bioassay example.

#### Solution
問題については教科書p.74参照。まず、$x_i$は薬の用量、$n_i$は動物、$y_i$は結果である。それぞれのグループi内において$y_i$は二項分布に従い、
$$
y_i|\theta_i \sim Bin(n_i, \theta_i)
$$
であり、$\theta_i$は致死率を表す。ここで、$\theta_i$と$x_i$の関係を以下のように書く。
$$
logit(\theta) = \alpha + \beta x_i
$$
この時、尤度は
$$
p(y|\theta) = \{logit^{-1}(\alpha + \beta x_i)\}^{y_i}\{1-logit^{-1}(\alpha + \beta x_i)\}^{n_i - y_i}
$$
対数尤度は
\begin{align}
c + y_i\{logit^{-1}(\alpha + \beta x_i)\} + (n_i-y_i)\{1-logit^{-1}(\alpha + \beta x_i)\}
\end{align}
ここでcは定数項。まずinformation matrixを求めるために、2階微分をそれぞれ考える。
\begin{align}
\frac{d^2\mathrm{log}p(y|\theta)}{d\alpha^2} = -\frac{n_iexp(\alpha + \beta x_i))}{(1+exp(\alpha + \beta x_i))^2}\\
\frac{d^2\mathrm{log}p(y|\theta)}{d\alpha d\beta} = -\frac{n_ix_iexp(\alpha + \beta x_i))}{(1+exp(\alpha + \beta x_i))^2}\\
\frac{d^2\mathrm{log}p(y|\theta)}{d\beta^2} = -\frac{n_ix_i^2exp(\alpha + \beta x_i))}{(1+exp(\alpha + \beta x_i))^2}
\end{align}
従って、情報行列は
\begin{align}
I(\hat\theta) = \begin{pmatrix}
\sum \frac{n_iexp(\alpha + \beta x_i))}{(1+exp(\alpha + \beta x_i))^2} & \sum\frac{n_ix_iexp(\alpha + \beta x_i))}{(1+exp(\alpha + \beta x_i))^2}\\
\sum\frac{n_ix_iexp(\alpha + \beta x_i))}{(1+exp(\alpha + \beta x_i))^2} & \sum\frac{n_ix_i^2exp(\alpha + \beta x_i))}{(1+exp(\alpha + \beta x_i))^2} \\
\end{pmatrix}
\end{align}
また、分散の近似値は情報行列の逆行列を求めればよい。
すなわち、$I(\hat\theta) = \begin{pmatrix} a & c\\ c&b \end{pmatrix}$とすれば、$\alpha$の分散は
\begin{aligned}
\frac{b}{ab-c^2}
\end{aligned}
同様に$\beta$の分散は
\begin{aligned}
\frac{a}{ab-c^2}
\end{aligned}

## 5

### a
#### Question
Suppose x and y are independent normally distributed random variables, where x has mean 4 and standard deviation 1, and y has mean 3 and standard deviation 2. What are the mean and standard deviation of y/x? Compute this using simulation.

#### Solution

set.seed(1234)
x <- rnorm(10000, mean = 4, sd = 1)
y <- rnorm(10000, mean = 3, sd = 2)
mean(y/x)
sd(y/x)

よって、平均は0.8003443、標準偏差は0.6223355

### b
#### Question
Suppose x and y are independent random variables, where x has mean 4 and standard deviation 1, and y has mean 3 and standard deviation 2. What are the approximate mean and standard deviation of y/x? Determine this without using simulation.
#### Solution

(辻村試案: 平均計算)
$E(x)=4,V(x)=1$より、$p(x)$は$4-\sqrt{3}\leq x \leq 4+\sqrt{3}$の範囲で$p(x)=\frac{1}{2\sqrt{3}}$となる一様分布に従う。
変数変換について、以下のように変換できる。([参考](https://ocw.u-tokyo.ac.jp/lecture_files/engin_05/2/notes/ja/ishikawa2.pdf#page=21))
\begin{align}
P(x<X,X<x+\Delta x)=P(z<Z,Z<z+\Delta z)\\
&\Leftrightarrow \int_{x}^{x+\Delta x}f(x')dx'&=\int_{z}^{z+\Delta z}g(z')dz'\\
&\Leftrightarrow f(x)dx = g(z)dz\\
&\Leftrightarrow g(z) = f(x)\frac{dx}{dz}
\end{align}

$g(z)=f(x)\frac{dx}{dz}$を用いると、$E(z)=\int{zg(z)dz}=\int{zf(x)dx}$と変換できる。
$a=4-\sqrt{3}, b=4+\sqrt{3},Z=1/X$とおいたとき、$a \leq x \leq b$より、
\begin{align}
E(\frac{1}{x})&= \int_{a}^{b}{\frac{1}{x}\frac{1}{b-a}dx}\\
&= \frac{1}{b-a}(\log{b} - \log{a})\\
&= \frac{1}{b-a}\log{\frac{b}{a}}\\
&= \frac{1}{2\sqrt{3}}\log{\frac{4+\sqrt{3}}{4-\sqrt{3}}}\\
\end{align}

$1/x$と$y$が独立であると仮定すると、$E(y)=3$より、
\begin{align}
E(\frac{y}{x}) &= E(y)E(\frac{1}{x})\\
&= \frac{3}{2\sqrt{3}}\log{\frac{4+\sqrt{3}}{4-\sqrt{3}}}\\
\end{align}

# 平均値
sqrt(3)/2*log(4+sqrt(3)/4-sqrt(3))
>[1] 0.8604897

### c
#### Question
What assumptions are required for the approximation in (b) to be reasonable?
#### Solution

## 6
Statistical decision theory: a decision-theoretic approach to the estimation of an unknown parameter θ introduces the loss function L(θ, a) which, loosely speaking, gives the cost of deciding that the parameter has the value a, when it is in fact equal to θ. The estimate a can be chosen to minimize the posterior expected loss,
$E(L(a|y))= \int L(\theta,a)p(\theta|y)d\theta$
This optimal choice of a is called a Bayes estimate for the loss function L. Show that:

### a
#### Question
If $L(θ, a) = (θ − a)^{2}$ (squared error loss), then the posterior mean, E(θ|y), if it exists, is the unique Bayes estimate of θ.
#### Solution
$L(a|y)$の微分を取る。
\begin{align}
\frac{d}{da}E(L(a|y)) &= \frac{d}{da}\int(\theta - a)^{2}p(\theta|y)d\theta\\
&=　 \frac{d}{da} \int \theta^{2} p(\theta|y)d\theta+ \frac{d}{da} \int (-2a\theta) p(\theta|y)d\theta+ \frac{d}{da} \int a^{2} p(\theta|y)d\theta \\
\\
&=0 - 2 \int\theta p(\theta|y)d\theta + 2a \int p(\theta|y)d\theta\\
& = -2(E(\theta|y) - a)
\end{align}

したがって、$a = E(\theta|y)$のとき、この微分は０を取る。二階微分は正となるのは明確なので、$L(a|y)$の最小値は$a = E(\theta|y)$のときである。したがって、Bayes estimateは$E(\theta|y)$となる。

### b
#### Question
If L(θ, a) = |θ − a|, then any posterior median of θ is a Bayes estimate of θ.
#### Solution
$k_{0} = k_{1} = 0$としたときの、ｃの回答の結果を適応すれば良い。

### c
#### Question
If $k_{0}$ and $k_{1}$ are nonnegative numbers, not both zero, and
\begin{align}
L(θ,a)=\begin{cases}
{k_{0}(θ−a) \ (θ≥a)}\\
{k_{1}(a−θ) \ (θ<a)}
\end{cases}
\end{align}
then any $\frac{k_{0}}{k_{0}+k_{1}}$ quantile of the posterior distribution p(θ|y) is a Bayes estimate of θ.

#### Solution
$L(a|y)$の微分を取る。
\begin{align}
\frac{d}{da}E(L(a|y)) &= \frac{d}{da}　( \int^{a}_{-\infty}k_{1}(a-\theta)p(\theta|y)d\theta + \int^{\infty}_{a}k_{0}(\theta-a)p(\theta|y)d\theta)\\
&=\frac{d}{da}(k_{1}(\int^{a}_{-\infty}a p(\theta|y)d\theta - \int^{a}_{-\infty}\theta p(\theta|y)d\theta) + k_{0}(\int^{\infty}_{a} \theta p(\theta|y)d\theta - \int^{\infty}_{a} a p(\theta|y)d\theta)) \\
&=　k_{1} \int^{a}_{-\infty} p(\theta|y)d\theta - k_{0} \int^{\infty}_{a} p(\theta|y)d\theta\\
&=　k_{1} \int^{a}_{-\infty} p(\theta|y)d\theta - k_{0} (1 - \int^{a}_{-\infty} p(\theta|y)d\theta)\\
& = (k_{1}+k_{0})\int^{a}_{-\infty} p(\theta|y)d\theta -k_{0}
\end{align}

これは$p(\theta|y) = \frac{k_{0}}{k_{0}+k_{1}}$のとき、０を取る。よって$\frac{k_{0}}{k_{0}+k_{1}}$はBayes estimateである。

## 4.7.
### Question
Unbiasedness: prove that the Bayesian posterior mean, based on a proper prior distri- bution, cannot be an unbiased estimator except in degenerate problems (see Bickel and Blackwell, 1967, and Lehmann, 1983, p. 244).

### Solution
事後分布の平均を $m(y) = E(\theta|y)$として、$m(y)$を$\theta$の推定量とする。
このとき、$m(y)$が$\theta$の不偏推定量であるためには、$E(m(y)|\theta)=\theta$である必要がある。
ここで、$\theta m(y)$ の周辺期待値は以下の通り表現できる。
\begin{align}
E(\theta m(y))=E[E(\theta m(y)|\theta)]=E[\theta^2]
\end{align}

また、異なる式変形として、
\begin{align}
E(\theta m(y))=E[E(\theta m(y)|\theta)]=E[m(y)^2]
\end{align}

とも表現できる。
以上から、
\begin{align}
E[\theta^2]+E[m(y)^2]&=2E(\theta m(y))\\
E[(m(y)-\theta)^2]&=0
\end{align}

この等式が満たされるのは確率1で$m(y)=\theta$となる場合のみである。

## 4.8.
### Question
Regression to the mean: work through the details of the example of mother’s and daughter’s heights on page 94, illustrating with a sketch of the joint distribution and relevant conditional distributions.

### Solution
p.94の議論は、成人女性の身長$\theta$とその母親の身長$y$を推定するというもの。
仮定として、$\theta$と$y$がそれぞれ同じ平均160、同じ分散であり、相関係数が0.5の場合を考えている。

\begin{align}
\mu = \begin{pmatrix}
160 \\
160 \\
\end{pmatrix}, \Sigma = \begin{pmatrix}
\sigma^2 & 0.5\sigma^2 \\
0.5\sigma^2 & \sigma^2 \\
\end{pmatrix}
\end{align}

library(mvtnorm)
library(scatterplot3d)

sigma <- matrix(c(1,0.5,0.5,1), ncol=2) # 分散共分散行列
x <- rmvnorm(n=10000, mean=c(160,160), sigma=sigma) # 2変数正規分布から乱数を生成
scatterplot3d(x[,1], x[,2], # 同時分布を3次元に描画
              dmvnorm(x, mean=c(160,160),
                      sigma=sigma), highlight=TRUE)


library(tidyverse)
x1 <- as_tibble(x)
colnames(x1) <- c("a","b")
x1 %>% filter(a>159,a<161) %>% .$b %>% hist(freq=FALSE)
x1 %>% filter(a>170,a<171) %>% .$b %>% hist(freq=FALSE)


## 13
### Question
Objections to Bayesian inference: discuss the criticism, ‘Bayesianism assumes: (a) Either
a weak or uniform prior [distribution], in which case why bother?, (b) Or a strong prior
[distribution], in which case why collect new data?, \(c\) Or more realistically, something
in between, in which case Bayesianism always seems to duck the issue’ (Ehrenberg, 1986).
Feel free to use any of the examples covered so far to illustrate your points.

### Solution
(a) Bayesian inferenceは興味のある変数を直接推定する方法だから、weak priorを使おうが問題はない。3.5節（分散既知の多変数正規モデルで平均を調べる）や、3.7節（Bioassay experimentで、薬の致死効果や半数致死量を調べる）の例に鑑みれば、Baysian inferenceは単刀直入な推定を行なっていることが分かる。これ以上良い方法があるか？

(b) Strong priorを設定するのに新しいデータを集める必要はない。Chapter 2の出生性比率の例を見れば分かるように、先行研究の理論や分析結果から引っ張ってくればいい話だからである。新しいデータは、他の条件で適用できるかを見極めるなど、よりモデルの一般化を目指す場合に集めればいいだけである。2.4節の前置胎盤の例では、strong priorを形成するためにドイツにおける全ての出生に関するデータを集めることはできない。しかし、新しいデータを得た場合に女児出生率が変わる可能性は考慮しておかなければならない。

\(c\) 「Bayesianが問題を避けている（duck the issue）」ーここでいう問題とは、priorとデータにはどのようなウェイトをかけるべきか、という問題である。Bayesianでは、データの利用可能性に基づきウェイトを設定している。Chapter 5で取り上げられる例では、このウェイトについての論点が議論される。そこでは、priorとdataのウェイトはdataそれ自身によって決定されていることが示される。Bayesianのどこが問題を避けているのか？

## 14
### Question
Objectivity and subjectivity: discuss the statement, ‘People tend to believe results that support their preconceptions and disbelieve results that surprise them. Bayesian methods encourage this undisciplined mode of thinking.’

### Solution
この命題には2つの問題がある。第一に、BDAが絶対的に主観的な結果の解釈を促すかという問題であり、第二に、BDAがその他の分析手法と比べて相対的に主観的な結果の解釈を促すかという問題である。

第一の点から考える。BDAにおいて主観の問題として大きいのは、事前分布の設定である。まず$n$が十分大きい場合から考える。この場合、事前分布は事後分布に対して大きく影響を与えない。例えば、3.10.2においてブッシュへの投票を二項分布でモデリングした際は、事前分布を$Beta(1, 1)$にすることで、尤度は$Beta(294, 364)$に対して、事後分布は、$Beta(295, 365)$であり、尤度の影響の方が遥かに大きいことがわかる。また、4.2で考察したように、$n \rightarrow \infty$の正規近似においては、事前分布は無視することができた。よって、$n$が十分大きい場合は、事前分布の影響力を低くすることを通じて、データに基づいた解釈が可能になる。また、$n$が小さい場合も、2.9で行ったように、informative priorと合わせてweakly informative priorを別途設定した分析を行うことにより、より客観的な分析を行うことが可能である。

第二に、そのほかの分析手法との比較から考える。まず、最尤推定の場合、$n\rightarrow \infty$において正規近似する場合、$y|\mu, \sigma \sim N(\mu, \sigma^2)$として、
\begin{align}
l(\theta|y) &\propto \log \prod \exp(-\frac{1}{2\sigma^2}(y_i - \theta)^2) \\
&= \sum(-\frac{1}{2\sigma^2}(y_i - \theta)^2)
\end{align}
よって、
\begin{align}
&\frac{d}{d\theta}l(\theta|y) \propto \frac{d}{d\theta} \sum(-\frac{1}{2\sigma^2}(y_i - \theta)^2) = 0 \\
&\Leftrightarrow \hat{\theta} = \bar{y}　&(\because \bar{y} = \sum y_i/n )
\end{align}
この値は、BDAにおいて、事前分布を一様分布にした場合に等しい。(一様でない場合は、p. 42)。よって、$n\rightarrow \infty$の場合は、BDAか最尤推定かは大きな問題ではない。なお、この近似ができない場合は、そもそも$\hat{\theta}$が十分統計量では無いので、最尤推定は非効率である。

次に、仮説検定との違いを考える。(著者の立場では)「$\theta =0$か?」という仮説よりも、$\theta$の事後分布を考えた方がより便利であろう。また、単にp-value < 0.05かどうかを見るよりも、「$\theta > 0$の事後確率が84%である」とした方がより情報に富んだ考察が可能である。

## 15
Coverage of posterior intervals:
a. Consider a model with scalar parameter $\theta$. Prove that, if you draw $\theta$ from the prior, draw $y|\theta$ from the data model, then perform Bayesian inference for $\theta$ given $y$, that there is a 50% probability that your 50% interval for $\theta$ contains the true value.
b. Suppose $\theta \sim N(0, 2^2)$ and $y|\theta \sim N(\theta, 1)$. Suppose the true value of $\theta$ is 1. What is the coverage of the posterior 50% interval for $\theta$? (You have to work this one out; it’s not 50% or any other number you could just guess.)
c. Suppose $\theta \sim N(0, 2^2)$ and $y|\theta \sim N(\theta, 1)$. Suppose the true value of $\theta$ is $\theta_0$. Make a plot showing the coverage of the posterior 50% interval for $\theta$, as a function of $\theta_0$.

### a.
$f(\theta|y) \propto f(y|\theta)f(\theta)$だから?

### b. & c.
\begin{align}
f(\theta|y) &\propto f(y|\theta)f(\theta) \\
&\propto \exp\left(-\frac{1}{2} (y - \theta)^2\right) \exp\left(-\frac{1}{2} \cdot \frac{\theta^2}{4}\right) \\
&\propto \exp\left(-\frac{1}{2} \cdot \frac{5}{4} \left(\theta^2 - 2 \cdot \frac{4}{5}y\theta\right)\right) \\
&\propto \exp\left(-\frac{1}{2} \cdot \frac{5}{4} \left(\theta - \frac{4}{5}y\right)^2\right) \\
&\propto N\left(\frac{4}{5}y, \frac{4}{5}\right).
\end{align}

以下は上式を用いて、$\theta$の事後確信区間が$\theta_0$を確率をシミュレートするものである。

bda3_4_7_15 <- function(theta0, seed = 123){
  set.seed(seed)
  # First, generate y
  theta <- rnorm(5000, 0, 2)
  y <- rnorm(5000, theta, 1)
  # Then, evaluate whether each posterior CI contains the theta0
  CIs <- matrix(NA, nrow = 5000, ncol = 2)
  colnames(CIs) <- c("lower", "upper")
  CIs[, 1] <- qnorm(0.25, (4*y)/5, sqrt(4/5))
  CIs[, 2] <- qnorm(0.75, (4*y)/5, sqrt(4/5))
  CIs <- as.data.frame(CIs)
  coverage <- ifelse(CIs$lower <= theta0 & CIs$upper >= theta0, 1, 0)
  sum(coverage) / 5000
}