BDA勉強会
      • Sharing URL Link copied
      • /edit
      • View mode
        • Edit mode
        • View mode
        • Book mode
        • Slide mode
        Edit mode View mode Book mode Slide mode
      • Customize slides
      • Note Permission
      • Read
        • Owners
        • Signed-in users
        • Everyone
        Owners Signed-in users Everyone
      • Write
        • Owners
        • Signed-in users
        • Everyone
        Owners Signed-in users Everyone
      • Engagement control Commenting, Suggest edit, Emoji Reply
    • Invite by email
      Invitee

      This note has no invitees

    • Publish Note

      Share your work with the world Congratulations! 🎉 Your note is out in the world Publish Note

      Your note will be visible on your profile and discoverable by anyone.
      Your note is now live.
      This note is visible on your profile and discoverable online.
      Everyone on the web can find and read all notes of this public team.
      See published notes
      Unpublish note
      Please check the box to agree to the Community Guidelines.
      View profile
    • Commenting
      Permission
      Disabled Forbidden Owners Signed-in users Everyone
    • Enable
    • Permission
      • Forbidden
      • Owners
      • Signed-in users
      • Everyone
    • Suggest edit
      Permission
      Disabled Forbidden Owners Signed-in users Everyone
    • Enable
    • Permission
      • Forbidden
      • Owners
      • Signed-in users
    • Emoji Reply
    • Enable
    • Versions and GitHub Sync
    • Note settings
    • Note Insights New
    • Engagement control
    • Make a copy
    • Transfer ownership
    • Delete this note
    • Insert from template
    • Import from
      • Dropbox
      • Google Drive
      • Gist
      • Clipboard
    • Export to
      • Dropbox
      • Google Drive
      • Gist
    • Download
      • Markdown
      • HTML
      • Raw HTML
Menu Note settings Note Insights Versions and GitHub Sync Sharing URL Help
Menu
Options
Engagement control Make a copy Transfer ownership Delete this note
Import from
Dropbox Google Drive Gist Clipboard
Export to
Dropbox Google Drive Gist
Download
Markdown HTML Raw HTML
Back
Sharing URL Link copied
/edit
View mode
  • Edit mode
  • View mode
  • Book mode
  • Slide mode
Edit mode View mode Book mode Slide mode
Customize slides
Note Permission
Read
Owners
  • Owners
  • Signed-in users
  • Everyone
Owners Signed-in users Everyone
Write
Owners
  • Owners
  • Signed-in users
  • Everyone
Owners Signed-in users Everyone
Engagement control Commenting, Suggest edit, Emoji Reply
  • Invite by email
    Invitee

    This note has no invitees

  • Publish Note

    Share your work with the world Congratulations! 🎉 Your note is out in the world Publish Note

    Your note will be visible on your profile and discoverable by anyone.
    Your note is now live.
    This note is visible on your profile and discoverable online.
    Everyone on the web can find and read all notes of this public team.
    See published notes
    Unpublish note
    Please check the box to agree to the Community Guidelines.
    View profile
    Engagement control
    Commenting
    Permission
    Disabled Forbidden Owners Signed-in users Everyone
    Enable
    Permission
    • Forbidden
    • Owners
    • Signed-in users
    • Everyone
    Suggest edit
    Permission
    Disabled Forbidden Owners Signed-in users Everyone
    Enable
    Permission
    • Forbidden
    • Owners
    • Signed-in users
    Emoji Reply
    Enable
    Import from Dropbox Google Drive Gist Clipboard
       Owned this note    Owned this note      
    Published Linked with GitHub
    • Any changes
      Be notified of any changes
    • Mention me
      Be notified of mention me
    • Unsubscribe
    # 2章章末問題 ## 2.11.1 ### Question Posterior inference: suppose you have a Beta(4, 4) prior distribution on the probability $\theta$ that a coin will yield a ‘head’ when spun in a specified manner. The coin is independently spun ten times, and ‘heads’ appear fewer than 3 times. You are not told how many heads were seen, only that the number is less than 3. Calculate your exact posterior density (up to a proportionality constant) for $\theta$ and sketch it. ### Solution 与えられたデータと事前分布から事後分布を求めれば良い。事前分布は \begin{align} p(\theta) \propto \theta^3 (1-\theta)^3 \end{align} 一方、事前分布からデータが発生する確率は、 \begin{align} p(y|\theta) &= \binom{10}{0} (1-\theta)^{10} + \binom{10}{1}\theta^1 (1-\theta)^9 + \binom{10}{2}\theta^2 (1-\theta)^8 \\ &= (1-\theta)^{10} + 10 \theta^1 (1-\theta)^9 + 45\theta^2 (1-\theta)^8 \end{align} よって事後分布は、 \begin{align} p(\theta|y) &\propto p(\theta)p(y|\theta) \\ &\propto \theta^3 (1-\theta)^3 ((1-\theta)^{10} + 10 \theta^1 (1-\theta)^9 + 45\theta^2 (1-\theta)^8) \\ &= \theta^{y+3} (1-\theta)^{13-y}((1-\theta)^{10} + 10 \theta^1 (1-\theta)^9 + 45\theta^2 (1-\theta)^8) \\ &= \theta^{3} (1-\theta)^{13} + 10 \theta^{4} (1-\theta)^{12} + 45\theta^{5} (1-\theta)^{11} \end{align} なお、現実の確率としては総和が1になるように定数をかける必要がある。 ```{r gegag} theta <- seq(0,1,.01) dens <- theta^3*(1-theta)^13 + 10*theta^4*(1-theta)^12 + 45*theta^5*(1-theta)^11 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) ``` ## 2.11.2 ### Question Predictive distributions: consider two coins, C1 and C2, with the following characteristics: Pr(heads|C1)=0.6 and Pr(heads|C2)=0.4. Choose one of the coins at random and imagine spinning it repeatedly. Given that the first two spins from the chosen coin are tails, what is the expectation of the number of additional spins until a head shows up? ### Solution コインの表が出る確率を、$\theta$とおくと、表が出るまでの回数を$x$として、その期待値は、 \begin{align} E(x|\theta) = \theta +2\theta(1-\theta)+3\theta(1-\theta)^2... \end{align} よって、 \begin{align} (1-\theta)E(x|\theta) = \theta(1-\theta) +2\theta(1-\theta)^2+3\theta(1-\theta)^3... \end{align} 今、両式の差をとると、 \begin{align} \theta E(x|\theta) &= \theta +\theta(1-\theta)+\theta(1-\theta)^2... \end{align} よって、 \begin{align} (1-\theta)\theta E(x|\theta) &= \theta(1-\theta) +\theta(1-\theta)^2+\theta(1-\theta)^3... \end{align} 今、以上の2つの式の差をとると、 \begin{align} \theta^2E(x|\theta) &= \theta \end{align} よって、$E(x|\theta) = 1/\theta$となる。 よって、表が出るまでの回数の期待値は、2回裏の観察結果を$y = TT$として、 \begin{align} E(x|y = TT) &= E(E(x|y = TT, \theta)|y=TT) \\ &= E(x|y = TT, \theta = 0.6)Pr(\theta = 0.6 |y=TT) + E(x|y = TT, \theta= 0.4)Pr(\theta = 0.4|y = TT) \\ &= 4/13 \cdot 1/0.6 + 9/13 \cdot 1/0.4 \\ &= 2.24 \end{align} ## 2.11.3 ### Question Predictive distributions: let y be the number of 6’s in 1000 rolls of a fair die. (a) Sketch the approximate distribution of y, based on the normal approximation. (b) Using the normal distribution table, give approximate 5%, 25%, 50%, 75%, and 95% points for the distribution of y. ### Solution 6の出る回数は二項分布に従うので、平均は \begin{align} E(y) &= 1000\cdot 1/6 \ &= 166.6 \end{align} そして、分散は、 \begin{align} V(y) &= 1000\cdot 1/6 \cdot 5/6\ &= 138.9 \end{align} よって、 ```{r gwegweg} y <- seq(120,220,.5) dens <- dnorm (y, 166.6, sqrt(138.9)) plot (y, dens, ylim=c(0,1.1*max(dens)), type="l", xlab="y", ylab="", xaxs="i", yaxs="i", yaxt="n", bty="n", cex=2) ``` 5%点は、 $166.7 − 1.65\cdot 11.8 = 147.2$ 25%点は、$166.7 − 0.67 \cdot 11.8 = 158.8$ 50%点は、$166.7$ 75%点は、$166.7 + 0.67 \cdot 11.8 = 174.6$ 95%点は、$166.7 + 1.65 \cdot 11.8 = 186.1$ 今、yは整数なので、147, 159, 167, 175, 186になる。 ## 2.11.4 ```{r gwegweg} # (a) density <- function(y, theta){ dnorm(y, 1000*theta, sqrt(1000*theta*(1-theta))) } y <- c(50:300) normal_approx <- 0.25*density(y, 1/12) + 0.5*density(y, 1/6) + 0.25*density(y, 1/4) plot(y, normal_approx, type = "l", ylab = "density") # (b) qnorm(0.2, 1000*(1/12), sqrt(1000*(1/12)*(11/12))) qnorm(0.9999, 1000*(1/12), sqrt(1000*(1/12)*(11/12))) qnorm(0.5, 1000*(1/6), sqrt(1000*(1/6)*(5/6))) qnorm(0.9999, 1000*(1/6), sqrt(1000*(1/6)*(5/6))) qnorm(0.8, 1000*(1/4), sqrt(1000*(1/4)*(3/4))) ``` ## 2.11.5 ### a. \begin{align} \Pr(y = k) &= \int^1_0 \Pr(y = k|\theta) d\theta \\ &= \binom{n}{k} \int^1_0 \theta^k (1 - \theta)^{n - k} d\theta \\ &= \binom{n}{k} \frac{\Gamma(k + 1)\Gamma(n - k + 1)}{\Gamma(n + 2)} \\ &= \frac{n!}{k!(n - k)!} \frac{k!(n - k)!}{(n + 1)!} \\ &= \frac{1}{n + 1} \end{align} ### b. 事後分布が$\mathrm{Beta}(y + \alpha, n - y + \beta)$なので、事後平均は$\frac{y + \alpha}{n + \alpha + \beta}$。また、事後平均が事前平均のデータの平均の間にあることから、0以上1以下のパラメーター$\lambda$を用いて、 \begin{align} \frac{y + \alpha}{n + \alpha + \beta} &= \lambda \frac{\alpha}{\alpha + \beta} + (1 - \lambda) \frac{y}{n} \\ &= \frac{y}{n} + \lambda (\frac{\alpha}{\alpha + \beta} - \frac{y}{n}) \\ \end{align} と書けるはずである。式を変形すると、 \begin{align} \frac{y + \alpha}{n + \alpha + \beta} - \frac{y}{n} &= \lambda (\frac{\alpha}{\alpha + \beta} - \frac{y}{n}) \\ \frac{n\alpha - y\alpha - n\beta}{(n + \alpha + \beta)n} &= \lambda (\frac{n\alpha - y\alpha - n\beta}{(\alpha + \beta)n}) \\ \lambda &= \frac{\alpha + \beta}{n + \alpha + \beta}. \\ \end{align} $\alpha >0$、$\beta > 0$、$n > 0$より、$0 < \lambda < 1$。 ### c. 事前分布の分散は$\frac{1}{12}$。事後分布は$\mathrm{Beta}(y + 1, n - y + 1)$なので、分散は \begin{align} \mathrm{var}(\theta|y) &= \frac{(y + 1)(n - y + 1)}{(n + 2)^2 (n + 3)} \\ &= \frac{(y + 1)}{(n + 2)} \frac{(n - y + 1)}{(n + 2)} \frac{1}{(n + 3)} \\ &< \frac{1}{4} \times \frac{1}{3} \end{align} となる。 ### d. $\alpha = 1$、$\beta = 10$、$n = 1$、$y = 1$と設定すると、事前分散は \begin{align} \frac{10}{11 \times 11 \times 12} \approx 0.0069, \end{align} 事後分散は \begin{align} \frac{20}{12 \times 12 \times 13} \approx 0.0107. \end{align} ## 2.11.6 \begin{align} y_j | \theta_j &\sim \mathrm{Poisson} (10n_j \theta_j) \\ \theta_j &\sim \mathrm{Gamma} (\alpha, \beta) \end{align} なので、式 (1.8) より \begin{align} \mathrm{E} (y_j) &= \mathrm{E} (\mathrm(E(y_j|\theta_j))) \\ &= \mathrm{E} (10n_j \theta_j) \\ &= 10n_j \mathrm{E}(\theta_j) \\ &= 10n_j \frac{\alpha}{\beta}. \end{align} また、式 (1.9) より \begin{align} \mathrm{var} (y_j) &= \mathrm{E} (\mathrm{var} (y_j|\theta_j)) + \mathrm{var} (\mathrm{E} (y_j|\theta_j)) \\ &= \mathrm{E} (10n_j \theta_j) + \mathrm{var} (10n_j \theta_j) \\ &= 10n_j \mathrm{E} (\theta_j) + (10n_j)^2 \mathrm{var} (\theta_J) \\ &= 10n_j \frac{\alpha}{\beta} + (10n_j)^2 \frac{\alpha}{\beta^2}. \end{align} ## 2.11.7 ### a. #### Question For the binomial likelihood, $y \sim Bin(n, \theta)$, show that $p(\theta) \propto \theta^{−1}(1 − \theta)^{−1}$ is the uniform prior distribution for the natural parameter of the exponential family. #### Solution 教科書p.36にあるように変形していくと、 \begin{align} y | \theta \sim \binom{n}{y}\theta^y (1 − \theta)^{n-y} &= \binom{n}{y}(1-\theta)^n (\frac{\theta}{1-\theta})^y\\ &= \binom{n}{y}(1-\theta)^n \mathrm{exp}(\mathrm{log}(\frac{\theta}{1-\theta})y) \end{align} とでき、自然パラメータがロジット変換の形になる。今、$p(\theta) \propto \theta^{−1}(1 − \theta)^{−1}$が自然パラメータ$\phi(\theta)$に対して一様な事前分布であることを示すので、$\phi$に一様分布を仮定する。すなわち、 \begin{align} p(\phi) = c \end{align} ここで、Jeffrey's invariance principleを用いて確率変数変換を行う。Jeffrey's invariance principleは母数の一対一変換(例えば$\phi =h(\theta)$)の不変性から無情報を定義していた(p.52参照)。すなわち、 \begin{align} p(\phi) = p(\theta)|\frac{d\theta}{d\phi}| \end{align} より、$\theta$を一様にするためには \begin{align} p(\theta) &= p(\phi)|\frac{d\phi}{d\theta}|\\ &= c \frac{d}{d\theta}\mathrm{log}(\frac{\theta}{1-\theta}) \\ &\propto \theta^{−1}(1 − \theta)^{−1} \end{align} ### b. #### Question Show that if $y = 0$ or $y = n$, the resulting posterior distribution is improper #### Solution (a)の結果より事前分布を \begin{align} p(\theta) = \theta^{−1}(1 − \theta)^{−1} = Beta(0,0) \end{align} とおく。このとき事後分布は \begin{align} p(\theta|y) &\propto \theta^{−1}(1 − \theta)^{−1} \binom{n}{y}\theta^y (1 − \theta)^{n-y}\\ &= \binom{n}{y}\theta^{y-1} (1 − \theta)^{n-y-1}\\ &\propto Beta(y, n-y) \end{align} となる。y = 0のとき事後分布は \begin{align} p(\theta | y = 0) \propto \theta^{-1} (1 − \theta)^{n-1} \end{align} また、$y=n$の時は \begin{align} p(\theta | y = n) \propto \theta^{n-1} (1 − \theta)^{-1} \end{align} であるので、$\theta$と$1-\theta$をいれかえれば$y=0$と$y=n$は等しくなることから、片方が積分しても1にならないことを示せば十分である。以下では$y=0$のケースを数学的帰納法を用いて証明する。 n = 0の時(つまり事前分布)、 \begin{align} \int_{0}^{1}p(\theta | y = 0, n = 0)d\theta &= \int_{0}^{1}\theta^{-1} (1 − \theta)^{-1}d\theta\\ &\geq \int_{0}^{1}\frac{1}{\theta}d\theta \\ &= [log \theta]_{0}^{1} = \infty \end{align} より成立。 次に$n = k-1$で積分した値が∞になっていたとする。すなわち、 \begin{align} \int_{0}^{1}p(\theta | y = 0, n = k-1)d\theta = \int_{0}^{1}\theta^{-1} (1 − \theta)^{k-2}d\theta = \infty \end{align} と仮定する。この時、$n=k$において、 \begin{align} \int_{0}^{1}p(\theta | y = 0, n = k)d\theta &= \int_{0}^{1}\theta^{-1} (1 − \theta)^{k-1}d\theta\\ &= \int_{0}^{1}\theta^{-1} (1 − \theta)^{k-2}(1 − \theta)d\theta\\ &= [(1 - \theta)^{k-2}(1 - \frac{\theta}{2})]_{0}^{1}+\int_0^1 \frac{(1 - \theta)^{k - 2}}{\theta}+ (k-2)(1 - \theta)^{k-3} - \frac{(1 - \theta)^{k-2}}{2}- (n - 2)\theta\frac{(1 - \theta)^{k-3}}{2} d\theta \\ &= c + \int_0^1 \frac{(1-\theta)^{k - 2}}{\theta}d\theta \end{align} より、cは$\infty$より小さい定数であることから、n=kにおいても積分すると$\infty$になる。従って、任意のnについて$p(\theta | y = 0)$はimproper distributionである。 ## 2.11.8 A random sample of n students is drawn from a large population, and their weights are measured. The average weight of the $n$ sampled students is $\bar{y} = 150$ pounds. Assume the weights in the population are normally distributed with unknown mean $\theta$ and known standard deviation 20 pounds. Suppose your prior distribution for $\theta$ is normal with mean 180 and standard deviation 40. ### a. #### Question Give your posterior distribution for $\theta$ (Your answer will be a function of n) #### Solution 事前分布及び尤度が正規分布なので、事後分布も正規分布となる。平均と分散は以下のようになる。 \begin{align} \mu = \mathbb E(\theta \mid \bar y) &= \frac{\frac{180}{1600} + \frac{150n}{400}}{\frac{1}{1600} + \frac{n}{400}} \\ &= \frac{60(3 + 10n)}{1600} \cdot \frac{1600}{1 + 4n} \\ &= \frac{60(3 + 10n)}{1 + 4n} \\ \frac{1}{\sigma^2} &= \frac{1}{1600} + \frac{n}{400} \\ &= \frac{1 + 4n}{1600} \end{align} ### b. #### Question A new student is sampled at random from the same population and has a weight of $\tilde{y}$ pounds. Give a posterior predictive distribution for $\tilde{y}$ (Your answer will still be a function of n). #### Solution データ$\tilde{y}\sim N(\theta, \sigma^2)$を追加した場合の事後予測分布は一般に、 \begin{align} p(\tilde{y}) &= \int_{-\infty}^{\infty}{p(\tilde{y}|\theta)p(\theta|y)d\theta}\\ &\propto \mathrm{exp}\frac{-\tilde{y}^2}{2\sigma^2} \int_{-\infty}^{\infty}{\mathrm{exp}({-\frac{\theta^2}{2}(\frac{1}{\sigma^2}+\frac{1}{\tau_1^2})})\mathrm{exp}(\theta(\frac{\tilde{y}}{\sigma^2} + \frac{\mu_1}{\tau_1^2}))}d\theta\\ &\sim \mathrm{exp}(\frac{\tilde{y}-2\mu_1 \tilde{y} - \frac{\mu_1^2\sigma^2}{\tau_1^2}}{2(\sigma^2+\tau_1^2)}) \end{align} となり、これも正規分布$N(\mu_1, \sigma^2+\tau_1^2)$に従う。 よって、 \begin{align} p(\tilde{y}|y) &= \int_{-\infty}^{\infty}{N(\tilde{y}|\mu, 20^2)N(\theta|\mu_n, \tau_n)}\\ &= N(\tilde{y}|\mu_n, 20^2+\tau_n^2)\\ &= N(\tilde{y}|\frac{60(3 + 10n)}{1 + 4n},\frac{1}{\frac{1 + 4n}{1600}}+20^2) \end{align} ### c. #### Question For n = 10, give a 95% posterior interval for $\theta$ and a 95% posterior predictive interval for $\tilde{y}$. #### Solution 95% posterior interval for $\theta$は \begin{align} 150.7 \pm 1.96(6.25) = [138, 163] \end{align} 95% posterior interval for $\tilde{y}$は \begin{align} 150.7 \pm 1.96(20.95) = [110, 192] \end{align} ### d. #### Question For n = 100, give a 95% posterior interval for $\theta$ and a 95% posterior predictive interval for $\tilde{y}$. #### Solution 以下のようにしてシミュレーションすることができる。 ```{r gwegweg} mu <- function(n) 60 * (3 + 10 * n) / (1 + 4 * n) sigma <- function(n) 40 / sqrt(1 + 4 * n) percentiles <- c(0.05, 0.95) theta_posterior_interval <- qnorm(percentiles, mu(100), sigma(10)) y_posterior_interval <- qnorm(percentiles, mu(100), sqrt(sigma(100)^2 + 400)) ``` この時、$\theta$の95%信頼区間は[140.5, 161]、$\tilde{y}$は[110.7, 189.5]になる。 ## 2.11.9 Setting parameters for a beta prior distribution: suppose your prior distribution for $\theta$, the proportion of Californians who support the death penalty, is beta with mean 0.6 and standard deviation 0.3. ### a. #### Question Determine the parameters ($\alpha$ and $\beta$) of your prior distribution. Sketch the prior density function. #### Solution ベータ分布の期待値が0.6,標準偏差が0.3なので \begin{align} \frac{\alpha}{\alpha+\beta}=0.6\\ \frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}=0.3\\ \end{align} これを解くと$\alpha = 1.0, \beta = \frac{2}{3}$ ### b. #### Question A random sample of 1000 Californians is taken, and 65% support the death penalty. What are your posterior mean and variance for $\theta$? Draw the posterior density function. #### Solution 1000人($=n$)のうち65%が賛成であることから、賛成の人数は650人($=y$)。反対は350人。より、事後分布の平均は \begin{align} \frac{\alpha + y}{\alpha+\beta+ y} = 0.65 \end{align} 分散は \begin{align} \frac{(\alpha + y)(\beta+n-y)}{(\alpha+\beta+n)^2(\alpha+\beta+n+1)} = 0.00023 \end{align} ## 2.11.10 Discrete sample spaces: suppose there are N cable cars in San Francisco, numbered sequentially from 1 to N. You see a cable car at random; it is numbered 203. You wish to estimate N ### a Assume geometric prior with mean 100: \begin{align} (1/100)(99/100)^{N-1} \end{align} 203号の列車を見たので N < 203は不可能になります。 \begin{cases} \frac{1}{N}&\text{for $N \geq 203$}\\ 0&\text{otherwise}\\ \end{cases} \begin{align} p(N|data) \propto \frac{1}{N}(1/100)(99/100)^{N-1} \text{where $N \geq 203$} \end{align} \begin{align} \propto \frac{1}{N}(99/100)^{N} \text{where $N \geq 203$} \end{align} ### b What are the posterior mean and standard deviation of N? (Sum the infinite series analytically or approximate them on the computer.) To calculate the maen and SD we need to calculate a normalizing constant $c$ \begin{align} p(N|data) = c\frac{1}{N}(99/100)^{N} \end{align} p(N|data) sums to 1 so \begin{align} \frac{1}{c} = \sum_{203}^{N}\frac{1}{N}(99/100)^{N} \end{align} ```{r} calculate_sum <-function(n){ (1/n)*((0.99)^n) } c_inv <- sum(calculate_sum(203:3000)) c <- c_inv^-1 ``` from this $c = 21.46829$. 平均値を計算すると \begin{align} p(N|data) = c\sum_{203}^{N}\frac{1}{N}(99/100)^{N} \\ = \frac{21.47 (0.99)^{203}}{1 − 0.99} \\ = 279.1 \end{align} ## 2.11.11 Computing with a nonconjugate single-parameter model: suppose y1,...,y5 are independent samples from a Cauchy distribution with unknown center $\theta$ and known scale 1: p(yi|θ) ∝ 1/(1 + (yi − θ)2). Assume, for simplicity, that the prior distribution for θ is uniform on [0, 100]. Given the observations (y1,...,y5) = (43, 44, 45, 46.5, 47.5): ### a Compute the unnormalized posterior density function, p(θ)p(y|θ), on a grid of points $θ = 0, \frac{1}{m} , \frac{2}{m} ,..., 100$ for some large integer m. Using the grid approximation, compute and plot the normalized posterior density function, p(θ|y), as a function of θ ```{r} dens <- function (y, th){ density <- c() for (i in 1:length(th)){ cauchy_draws <- dcauchy(y, th[i], 1) cauchy_products <- prod(cauchy_draws) density <- c(density, cauchy_products) } return(density) } y <- c(43, 44, 45, 46.5, 47.5) theta <- seq(1, 100, 0.5) density <- dens(y,theta) density_normalized <- density/(0.5*sum(density)) ``` ![](https://i.imgur.com/YXLETR8.png) ### b Sample 1000 draws of θ from the posterior density and plot a histogram of the draws. ```{r} thetas <- sample(theta, 1000, density_normalized, replace=TRUE) hist(thetas) ``` ![](https://i.imgur.com/EHfvB4m.png) ### c Use the 1000 samples of θ to obtain 1000 samples from the predictive distribution of a future observation, y6, and plot a histogram of the predictive draws ```{r} new_obs <- rcauchy(length(thetas), thetas, 1) hist(new_obs) ``` ![](https://i.imgur.com/8iI41hr.png) ## 2.11.12 Jeffreys’ prior distributions: suppose y|θ ∼ Poisson(θ). Find Jeffreys’ prior density for θ, posson pdf is : \begin{align} p(y|\theta) = \frac{\theta^{y}e^{-\theta}}{y!} \end{align} Jeffrey's prior is \begin{align} J(\theta) = E(\frac{-d^{2}\text{log }p(y|\theta)}{d\theta^{2}}|\theta) \\ = E(y/\theta^{2}) \\ = 1/\theta \end{align} ## 2.11.13 Worldwide airplance fatalities | Year | Fatal accidents | Passenger death |Death rate (passenger deaths per 100 million passenger miles) | | -------- | -------- | -------- |------| | 1976 | 24 | 734 |0.19 | | 1977 | 25 | 516 |0.12 | | 1978 | 31 | 754 |0.15 | | 1979 | 31 | 877 |0.16 | | 1980 | 22 | 814 |0.14 | | 1981 | 21 | 362 |0.06 | | 1982 | 26 | 764 |0.13 | | 1983 | 20 | 809 |0.13 | | 1984 | 16 | 223 |0.03 | | 1985 | 22 | 1066 |0.15 | ### (a) #### Question Assume that the numbers of fatal accidents in each year are independent with a Poisson(θ) distribution. Set a prior distribution for θ and determine the posterior distribution based on the data from 1976 through 1985. Under this model, give a 95% predictive interval for the number of fatal accidents in 1986. You can use the normal approximation to the gamma and Poisson or compute using simulation. #### Solution $y_{i}$をyear *i* におけるfatal accidentの数、$\theta$をある1年間での事故数の期待値とする。 \begin{align} y_{i}| \theta \sim Poisson (\theta) \end{align} と表される。 ポワソン分布の共役事前分布はガンマ分布。前回の発表より、仮に$Gammma(\alpha, \beta)$を事前分布として用いたとき、事後分布は$Gammma(\alpha+\sum_{i=1}^{n} y_{i}, \beta + n)$となる。 事前分布を定める。ここでどうハイパーパラメーター$\alpha$、$\beta$を定めれば良いかは正直わからなかったが、解答では無情報事前分布を用いており、そのパラメーターは$\alpha=0$、$\beta=0$であった。$n=10$とここでは十分大きいので、無情報事前分布を用いることは正当化されるらしい。 よって、$Gamma(0,0)$を事前分布として用いる。fatal accidentの10年間での合計は238のため、事後分布は$Gamma(238, 10)$となる。 1986年における未知だが観測可能な事故発生数を$\tilde{y}$とおくと、 \begin{align} \tilde{y} \sim Poisson(\theta) \end{align} と表すことができ、ここでのパラメーター$\theta$は先ほどのガンマ分布$Gamma(238, 10)$に従う。 Rで以下のようにシミュレーションをすることで、$\tilde{y}$の95%予測信用区間は求めることができ、[14, 34]となる。 ```{r gwegweg} set.seed(1234) theta <- rgamma(1000, 238, 10) y1986 <- rpois(1000, theta) sort(y1986)[c(25, 976)] ``` ### (b) #### Question Assume that the numbers of fatal accidents in each year follow independent Poisson distributions with a constant rate and an exposure in each year proportional to the number of passenger miles flown. Set a prior distribution for θ and determine the posterior distribution based on the data for 1976–1985. (Estimate the number of passenger miles flown in each year by dividing the appropriate columns of Table 2.2 and ignoring round-off errors.) Give a 95% predictive interval for the number of fatal accidents in 1986 under the assumption that 8 × 1011 passenger miles are flown that year. #### Solution (Passenger Death/Death Rate)x(100 million miles)で各年のPassenger Miles (各旅客が飛行した距離の合計) が出せる。 | Year | Estimated number of passenger miles | | -------- | -------- | | 1976 | $3.863 \times 10^{11}$ | | 1977 | $4.300 \times 10^{11}$ | | 1978 | $5.027 \times 10^{11}$ | | 1979 | $5.481 \times 10^{11}$ | | 1980 | $5.814 \times 10^{11}$ | | 1981 | $6.033 \times 10^{11}$ | | 1982 | $5.877 \times 10^{11}$ | | 1983 | $6.223 \times 10^{11}$ | | 1984 | $7.433 \times 10^{11}$ | | 1985 | $7.106 \times 10^{11}$ | 前回のに説明変数を付け足したバージョンを考える $y_{i}$をyear *i* におけるfatal accidentの数、$\theta$をある1年間での事故数の期待値、$x_{i}$をyear *i* における旅客マイル(Passenger Mile)とする。このとき \begin{align} y_{i}| \theta \sim Poisson (x_{i}\theta) \end{align} と表される。 ポワソン分布の共役事前分布はガンマ分布。前回の発表より、仮に$Gammma(\alpha, \beta)$を事前分布として用いたとき、事後分布は$Gammma(\alpha+\sum_{i=1}^{n} y_{i}, \beta+\sum_{i=1}^{n} x_{i})$となる。 事前分布を定める。前回同様$Gamma(0,0)$を事前分布として用いると、事後分布は \begin{align} \theta|y \sim Gamma(238, 5.716 \times 10^{12}) \end{align} となる。 1986年における未知だが観測可能な事故発生数を$\tilde{y}$とおくと、 \begin{align} \tilde{y} \sim Poisson(8\times 10^{12}\theta) \end{align} と表すことができ、ここでのパラメーター$\theta$は先ほどのガンマ分布$Gamma(238, 5.716 \times 10^{12})$に従う。 Rで以下のようにシミュレーションをすることで、$\tilde{y}$の95%予測信用区間は求めることができ、[22, 46]となる。 ```{r gwegweg} set.seed(1234) theta <- rgamma(1000, 238, 5.716e12) y1986 <- rpois(1000, 8e11*theta) sort(y1986)[c(25, 976)] ``` ### (c ) #### Question Repeat (a) above, replacing ‘fatal accidents’ with ‘passenger deaths.’ #### Solution (a)と同様のことを"passenger death"で行う。 $y_{i}$をyear *i* におけるpassenger deathsの数、$\theta$をある1年間での事故数の期待値とする。 \begin{align} y_{i}| \theta \sim Poisson (\theta) \end{align} と表される。 ポワソン分布の共役事前分布はガンマ分布。前回の発表より、仮に$Gammma(\alpha, \beta)$を事前分布として用いたとき、事後分布は$Gammma(\alpha+\sum_{i=1}^{n} y_{i}, \beta + n)$となる。 $Gamma(0,0)$を事前分布として用いる。passenger deathの合計は6919なので、事後分布は$Gamma(6919, 10)$となる。 1986年における未知だが観測可能な事故発生数を$\tilde{y}$とおくと、 \begin{align} \tilde{y} \sim Poisson(\theta) \end{align} と表すことができ、ここでのパラメーター$\theta$は先ほどのガンマ分布$Gamma(6919, 10)$に従う。 Rで以下のようにシミュレーションをすることで、$\tilde{y}$の95%予測信用区間は求めることができ、[639, 747]となる。 ```{r gwegweg} set.seed(1234) theta <- rgamma(1000, 6919, 10) y1986 <- rpois(1000, theta) sort(y1986)[c(25, 976)] ``` ### (d) #### Question Repeat (b) above, replacing ‘fatal accidents’ with ‘passenger deaths.’ #### Solution (b)と同様のことを‘passenger deaths’におきかえて行う。 $y_{i}$をyear *i* におけるpassenger deathsの数、$\theta$をある1年間での事故数の期待値、$x_{i}$をyear *i* における旅客マイル(Passenger Mile)とする。このとき \begin{align} y_{i}| \theta \sim Poisson (x_{i}\theta) \end{align} と表される。 ポワソン分布の共役事前分布はガンマ分布。仮に$Gammma(\alpha, \beta)$を事前分布として用いたとき、事後分布は$Gammma(\alpha+\sum_{i=1}^{n} y_{i}, \beta+\sum_{i=1}^{n} x_{i})$となる。 事前分布を定める。前回同様$Gamma(0,0)$を事前分布として用いると、事後分布は \begin{align} \theta|y \sim Gamma(6919, 5.716 \times 10^{12}) \end{align} となる。 1986年における未知だが観測可能な事故発生数を$\tilde{y}$とおくと、 \begin{align} \tilde{y} \sim Poisson(8\times 10^{12}\theta) \end{align} と表すことができ、ここでのパラメーター$\theta$は先ほどのガンマ分布$Gamma(6919, 5.716 \times 10^{12})$に従う。 Rで以下のようにシミュレーションをすることで、$\tilde{y}$の95%予測信用区間は求めることができ、[908, 1032]となる。 ```{r gwegweg} set.seed(1234) theta <- rgamma(1000, 6919, 5.716e12) y1986 <- rpois(1000, 8e11*theta) sort(y1986)[c(25, 976)] ``` ### (e) #### Question In which of the cases (a)–(d) above does the Poisson model seem more or less reasonable? Why? Discuss based on general principles, without specific reference to the numbers in Table 2.2. Incidentally, in 1986, there were 22 fatal accidents, 546 passenger deaths, and a death rate of 0.06 per 100 million miles flown. We return to this example in Exercises 3.12, 6.2, 6.3, and 8.14. #### Solution ## 2.11.14 ### (a) #### Question Fill in the steps to derive (2.9)–(2.10), and (2.11)–(2.12). #### Solution 教科書p40の(2.9) (2.10)とは、分散が既知の正規分布から生じる状況において、事前分布を$p(\theta) \propto exp(-\frac{1}{2\tau_{0}^{2}}(\theta - \mu_{0})^{2})$ (,where $\theta \sim N(\mu_{0}, \tau_{0})$)とした時の事後分布が \begin{align} p(\theta | y) \propto exp(-\frac{1}{2}(\frac{(y -\theta)^{2}}{\sigma^{2}} + \frac{(\theta - \mu_{0})^{2}}{\tau_{0}^{2}})) \end{align} となり、これが \begin{align} p(\theta | y) \propto exp(-\frac{1}{2\tau_{1}^{2}}(\theta - \mu_{1})^{2}) \end{align} that is $\theta|y \sim N(\mu_{1}, \tau_{1}^{2})$ ,where \begin{align} \mu_{1} = \frac{\frac{1}{\tau_{0}^{2}} \mu_{0} + \frac{1}{\sigma^{2}}y} {\frac{1}{\tau_{0}^{2}} + \frac{1}{\sigma^{2}}},\ \frac{1}{\tau_{1}^{2}} = \frac{1}{\tau_{0}^{2}} + \frac{1}{\sigma^{2}} \end{align} と変形できることを指しているが、前者から後者が導かれるまでのステップを示せば良い。 \begin{align} p(\theta | y) &\propto p(\theta)p(y | \theta)\\ &\propto exp(-\frac{1}{2}(\frac{(\theta - \mu_{0})^{2}}{\tau_{0}^{2}})) exp(-\frac{1}{2}(\frac{(y -\theta)^{2}}{\sigma^{2}})\\ &= exp(-\frac{1}{2}(\frac{(y -\theta)^{2}}{\sigma^{2}} + \frac{(\theta - \mu_{0})^{2}}{\tau_{0}^{2}}))\\ &= exp(-\frac{1}{2}(\frac{y^{2} - 2y + \theta^{2}}{\sigma^{2}} + \frac{\theta^{2} - 2\mu_{0}\theta + \mu_{0}^{2}}{\tau_{0}^{2}}))\\ &= exp(-\frac{1}{2} \{ \theta^{2} (\frac{1}{\sigma^{2}} + \frac{1}{\tau_{0}^{2}}) - 2 \theta (\frac{y}{\sigma^{2}} + \frac{\mu_{0}}{\tau_{0}^{2}}) + (\frac{y^{2}}{\sigma^{2}} + \frac{\mu_{0}^{2}}{\tau_{0}^{2}})\})\\ &\propto exp(-\frac{1}{2} \{\theta^{2} (\frac{1}{\sigma^{2}} + \frac{1}{\tau_{0}^{2}}) - 2 \theta (\frac{y}{\sigma^{2}} + \frac{\mu_{0}}{\tau_{0}^{2}})\})\\ &= exp(-\frac{1}{2} (\frac{1}{\sigma^{2}} + \frac{1}{\tau_{0}^{2}})( \theta^{2} - 2 \theta \frac{\frac{y}{\sigma^{2}} + \frac{\mu_{0}}{\tau_{0}^{2}}}{\frac{1}{\sigma^{2}} + \frac{1}{\tau_{0}^{2}}}))\\ &\propto exp(-\frac{1}{2} (\frac{1}{\sigma^{2}} + \frac{1}{\tau_{0}^{2}})( \theta - \frac{\frac{y}{\sigma^{2}} + \frac{\mu_{0}}{\tau_{0}^{2}}}{\frac{1}{\sigma^{2}} + \frac{1}{\tau_{0}^{2}}})^{2})\\ &= exp(-\frac{1}{2\tau_{1}^{2}}(\theta - \mu_{1})^{2}) \end{align} 教科書p41、42の(2.11)(2.12)は複数のデータについて上記と同様のことを行なっており \begin{align} p(\theta | \boldsymbol{y}) &\propto p(\theta)p(\boldsymbol{y} | \theta)\\ &= p(\theta)\prod_{i=1}^{n} p(y_{i} | \theta) \\ &\propto exp(-\frac{1}{2\tau_{0}^{2}}(y_{i} - \theta)^{2}) \prod_{i=1}^{n}exp(-\frac{1}{2\tau_{0}^{2}}(\theta - \mu_{0})^{2}) \\ &\propto exp(-\frac{1}{2}((\frac{(\theta - \mu_{0})^{2}}{\tau_{0}^{2}}) + \frac{1}{\sigma^{2}}\sum_{i=1}^{n}(y_{i} -\theta)^{2}))\\ &= exp(-\frac{1}{2}((\frac{\theta^{2} - 2\theta \mu_{0} + \mu_{0}^{2}}{\tau_{0}^{2}}) + \frac{1}{\sigma^{2}}(\sum_{i=1}^{n} y_{i}^{2} - 2 \theta n \bar{y} + n\theta^{2})))\\ &\propto exp(-\frac{1}{2}(\theta^{2}(\frac{1}{\tau_{0}^{2}} + \frac{n}{\sigma^{2}}) -2\theta (\frac{\mu_{0}}{\tau_{0}^{2}} + \frac{n\bar{y}}{\sigma^{2}})))\\ &= exp(-\frac{1}{2}(\frac{1}{\tau_{0}^{2}} + \frac{n}{\sigma^{2}})(\theta^{2} -2\theta \frac{\frac{\mu_{0}}{\tau_{0}^{2}} + \frac{n\bar{y}}{\sigma^{2}}}{\frac{1}{\tau_{0}^{2}} + \frac{n}{\sigma^{2}}}))\\ &\propto exp(-\frac{1}{2}(\frac{1}{\tau_{0}^{2}} + \frac{n}{\sigma^{2}})(\theta - \frac{\frac{\mu_{0}}{\tau_{0}^{2}} + \frac{n\bar{y}}{\sigma^{2}}}{\frac{1}{\tau_{0}^{2}} + \frac{n}{\sigma^{2}}})^{2})\\ &= exp(-\frac{1}{2 \tau_{n}^{2}}(\theta- \mu_{n})^{2}) \end{align} , where $\mu_{n} = \frac{\frac{\mu_{0}}{\tau_{0}^{2}} + \frac{n}{\sigma^{2}}\bar{y}}{\frac{1}{\tau_{0}^{2}} + \frac{n}{\sigma^{2}}} , \frac{1}{\tau_{n}^{2}} = \frac{1}{\tau_{0}^{2}} + \frac{n}{\sigma^{2}}$ ### (b) #### Question Derive (2.11) and (2.12) by starting with a $N(\mu_{0}, \tau_{0}^{2})$ prior distribution and adding data points one at a time, using the posterior distribution at each step as the prior distribution for the next. #### Solution (a)より一回データを追加した時の事後分布は $\theta|y \sim N(\mu_{1}, \tau_{1}^{2})$となる。 二回目データを追加するとき、$N(\mu_{1}, \tau_{1}^{2})$を事前分布としたときの事後分布は、 \begin{align} N(\mu_{2}, \tau_{2}^{2}) \end{align} , where \begin{align} \mu_{2} = \frac{\frac{1}{\tau_{1}^{2}} \mu_{1} + \frac{1}{\sigma^{2}}y_{2}} {\frac{1}{\tau_{1}^{2}} + \frac{1}{\sigma^{2}}} = \frac{\frac{1}{\tau_{0}^{2}} \mu_{0} + \frac{1}{\sigma^{2}}(y_{1}+y_{2})} {\frac{1}{\tau_{0}^{2}} + \frac{2}{\sigma^{2}}},\\ \frac{1}{\tau_{2}^{2}} = \frac{1}{\tau_{1}^{2}} + \frac{1}{\sigma^{2}} = \frac{1}{\tau_{0}^{2}} + \frac{2}{\sigma^{2}} \end{align} i回目のデータを追加した時、$\theta | \boldsymbol{y} \sim N(\mu_{i}, \tau_{i}^{2})$が成立していると仮定する。 i+1回目データを追加するとき、$N(\mu_{i}, \tau_{i}^{2})$を事前分布としたときの事後分布は、 \begin{align} N(\mu_{i+1}, \tau_{i+1}^{2}) \end{align} , where \begin{align} \mu_{i+1} &= \frac{\frac{1}{\tau_{i}^{2}} \mu_{i} + \frac{1}{\sigma^{2}}y_{i+1}} {\frac{1}{\tau_{i}^{2}} + \frac{1}{\sigma^{2}}} = \frac{\frac{1}{\tau_{0}^{2}} \mu_{0} + \frac{1}{\sigma^{2}}\sum_{k=1}^{i}y_{k} + \frac{1}{\sigma^{2}}(y_{i+1})} {\frac{1}{\tau_{i}^{2}} + \frac{i+1}{\sigma^{2}}} = \frac{\frac{1}{\tau_{0}^{2}} \mu_{0} + \frac{1}{\sigma^{2}}(i+1)\bar{y}} {\frac{1}{\tau_{i}^{2}} + \frac{i+1}{\sigma^{2}}},\\ \frac{1}{\tau_{i+1}^{2}} &= \frac{1}{\tau_{i}^{2}} + \frac{1}{\sigma^{2}} = \frac{1}{\tau_{0}^{2}} + \frac{i}{\sigma^{2}}+ \frac{1}{\sigma^{2}} = \frac{1}{\tau_{0}^{2}} + \frac{i+1}{\sigma^{2}} \end{align} 従って、$n=1, i+1$でも成立したため、一般項$n$でも \begin{align} \theta | \boldsymbol{y} \sim N(\mu_{n}, \tau_{n}^{2}) \end{align} , where $\mu_{n} = \frac{\frac{\mu_{0}}{\tau_{0}^{2}} + \frac{n}{\sigma^{2}}\bar{y}}{\frac{1}{\tau_{0}^{2}} + \frac{n}{\sigma^{2}}} , \frac{1}{\tau_{n}^{2}} = \frac{1}{\tau_{0}^{2}} + \frac{n}{\sigma^{2}}$ も成立し、 順にデータを足した時でも(2.11)(2.12)は成立する。 ## 2.11.15 ### Question Beta distribution: assume the result, from standard advanced calculus, that \begin{align} \int_{0}^{1} u^{\alpha - 1}(1-u)^{\beta -1}du = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha + \beta)} \end{align} If $Z$ has a beta distribution with parameters $\alpha$ and $\beta$, find $E[Z^{m}(1 − Z)^{n}]$ for any non-negative integers $m$ and $n$. Hence derive the mean and variance of $Z$. ### Solution Suppose $Z \sim beta(\alpha, \beta)$. \begin{align} E[Z^{m}(1 − Z)^{n}] &= \frac{1}{B(\alpha, \beta)} \int_{0}^{1} z^{m}(1-z)^{n}z^{\alpha - 1}(1-z)^{\beta - 1} dz\\ &= \frac{1}{B(\alpha, \beta)} \int_{0}^{1} z^{m + \alpha -1}(1-z)^{n + \beta -1} dz\\ &= \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \frac{\Gamma (m + \alpha) \Gamma (n + \beta)}{\Gamma (m + n + \alpha + \beta)} \end{align} (メモ$\frac{(\alpha + \beta -1)!}{(\alpha - 1)!(\beta - 1)!}\frac{(m + \alpha -1)! (n + \beta -1)!}{(m + n + \alpha + \beta - 1)!}$ $\frac{\int_{0}^{\infty} t^{\alpha + \beta -1}e^{-t}dt}{\int_{0}^{\infty} t^{\alpha -1}e^{-t}dt \int_{0}^{\infty} t^{\beta -1}e^{-t}dt} \frac{\int_{0}^{\infty} t^{m + \alpha -1}e^{-t}dt \int_{0}^{\infty} t^{n + \beta -1}e^{-t}dt}{\int_{0}^{\infty} t^{m + n + \alpha + \beta - 1}e^{-t}dt}$) $Z$の期待値は上記の式の$m=1, n=0$とした時。 \begin{align} E[Z] &= \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \frac{\Gamma (1 + \alpha) \Gamma (\beta)}{\Gamma (\alpha + \beta + 1)} = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)} \frac{\Gamma (1 + \alpha)}{\Gamma (\alpha + \beta + 1)} \end{align} ここで、ガンマ関数の漸化式($\Gamma(s + 1) = s\Gamma(s)$)を用いて \begin{align} & \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)} \frac{\Gamma (1 + \alpha)}{\Gamma (\alpha + \beta + 1)}\\ &= \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)} \frac{\alpha \Gamma (\alpha)}{(\alpha + \beta)\Gamma (\alpha + \beta)}\\ &= \frac{\alpha}{\alpha + \beta} \end{align} 次に分散を求める。$Z^{2}$の期待値は先ほどの式を$m=2, n=0$とした時。 \begin{align} E[Z^{2}] &= \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \frac{\Gamma (2 + \alpha) \Gamma (\beta)}{\Gamma (\alpha + \beta + 2)} = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)} \frac{\Gamma (2 + \alpha)}{\Gamma (\alpha + \beta + 2)} \end{align} 同様にガンマ関数の漸化式を繰り返し用いて \begin{align} & \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)} \frac{\Gamma (2 + \alpha)}{\Gamma (\alpha + \beta + 2)}\\ & = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)} \frac{(1 + \alpha) \Gamma (1 + \alpha)}{(\alpha + \beta + 1)\Gamma (\alpha + \beta + 1)}\\ & = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)} \frac{(1 + \alpha) \alpha \Gamma ( \alpha)}{(\alpha + \beta + 1)(\alpha + \beta) \Gamma (\alpha + \beta)}\\ & = \frac{(1 + \alpha) \alpha}{(\alpha + \beta + 1)(\alpha + \beta)}\\ \end{align} 従って、分散は \begin{align} V(Z) &= E[Z^{2}] - E[Z]^{2}\\ &= \frac{(1 + \alpha) \alpha}{(\alpha + \beta + 1)(\alpha + \beta)} - \frac{\alpha^{2}}{(\alpha + \beta)^{2}}\\ &= \frac{\alpha (1 + \alpha) (\alpha + \beta) - \alpha^{2}(\alpha + \beta +1)}{(\alpha + \beta + 1)(\alpha + \beta)^{2}}\\ &= \frac{\alpha\beta}{(\alpha + \beta + 1)(\alpha + \beta)^{2}} \end{align} ## 2.11.16 ### Question Beta-binomial distribution and Bayes’ prior distribution: suppose y has a binomial dis- tribution for given n and unknown parameter θ, where the prior distribution of θ is Beta(α, β). (a) Find p(y), the marginal distribution of y, for y = 0, . . . , n (unconditional on θ). This discrete distribution is known as the beta-binomial, for obvious reasons. ### Solution \begin{align} p(y) &= \int p(y|\theta)p(\theta)d\theta \\ &= \int_{0}^{1}\binom{n}{y} \theta^y(1-\theta)^{n-y}Beta(\alpha, \beta) d\theta \\ &= \int_{0}^{1}\binom{n}{y} \theta^y(1-\theta)^{n-y}\frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta^{\alpha-1}(1-\theta)^{\beta-1} d\theta \\ &= \frac{\Gamma(n+1)}{\Gamma(y+1)\Gamma(n-y+1)}\frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)}\frac{\Gamma(y+\alpha)\Gamma(n-y+\beta)}{\Gamma(n+\alpha+\beta)} \end{align} ### Question (b) Show that if the beta-binomial probability is constant in y, then the prior distribution has to have α = β = 1. ### Solution \begin{align} p(y) \propto \frac{\Gamma(y+\alpha)\Gamma(n-y+\beta)}{\Gamma(y+1)\Gamma(n-y+1)} \end{align} ここで、 $\alpha=\beta=1$であるならば \begin{align} \frac{\Gamma(y+\alpha)\Gamma(n-y+\beta)}{\Gamma(y+1)\Gamma(n-y+1)} &= \frac{\Gamma(y+1)\Gamma(n-y+1)}{\Gamma(y+1)\Gamma(n-y+1)}\\ &=1 \end{align} となり、$y$の値に関わりなく一定 また、$p(y)$が$y$の値に関わりなく一定である場合、 $p(0)=p(1), p(0)=p(n)$が成立する。 $p(0)=p(n)$について、 \begin{align} \frac{\Gamma(\alpha)\Gamma(n+\beta)}{\Gamma(1)\Gamma(n+1)} &= \frac{\Gamma(n+\alpha)\Gamma(\beta)}{\Gamma(n+1)\Gamma(1)}\\ {\Gamma(\alpha)\Gamma(n+\beta)}&= \Gamma(n+\alpha)\Gamma(\beta) \end{align} ここで、$\Gamma(t)=(t-1)\Gamma(t-1)$より、 \begin{align} \Gamma(\beta+n)&= (\beta+n-1)(\beta+n-2)...(\beta+1)\beta\Gamma(\beta) \end{align} したがって、 \begin{align} \Gamma(\alpha)\Gamma(n+\beta)&= \Gamma(n+\alpha)\Gamma(\beta)\\ \Gamma(\alpha)(\beta+n-1)(\beta+n-2)...(\beta+1)\beta\Gamma(\beta)&=\Gamma(\beta)(\alpha+n-1)(\alpha+n-2)...(\alpha+1)\alpha\Gamma(\alpha) \end{align} 全ての項は正なので、これを満たすのは$\alpha=\beta$のとき $p(0)=p(1)$について、同様に$\Gamma(t)=(t-1)\Gamma(t-1)$を利用して \begin{align} \frac{\Gamma(\alpha)\Gamma(n+\beta)}{\Gamma(1)\Gamma(n+1)} &= \frac{\Gamma(1+\alpha)\Gamma(n-1+\beta)}{\Gamma(2)\Gamma(n)}\\ \frac{\Gamma(\alpha)(n+\beta-1)\Gamma(n+\beta-1)}{\Gamma(1)n\Gamma(n)} &= \frac{\alpha\Gamma(\alpha)\Gamma(n+\beta-1)}{1\Gamma(1)\Gamma(n)}\\ n+\beta-1&=n\alpha \end{align} $\alpha=\beta$より、 \begin{align} \alpha+n-1=n\alpha \end{align} これを満たすのは$\alpha=1$のときのみであるから、$\alpha=\beta=1$$\alpha=\beta=1$ したがって、$\alpha=\beta=1$は$p(y)$が$y$について一定であることの必要十分条件となる。 ## 2.11.17 Posterior intervals: unlike the central posterior interval, the highest posterior interval is not invariant to transformation. For example, suppose that, given $σ^2$, the quantity $nv/σ^2$ is distributed as $χ^2_n$, and that σ has the (improper) noninformative prior density p(σ) ∝ σ−1,σ > 0. ### Question (a) Prove that the corresponding prior density for $σ^2$ is $p(σ^2) ∝ σ^{−2}$. ### Solution $\delta^2=u$とおく \begin{align} p(\delta^2)=p(u)&=\frac{p(\sqrt{u})d\sqrt{u}}{du}\\ &=p(\delta)(\frac{1}{2})u^{-2}\\ &=\frac{1}{2}p(\delta)\frac{1}{\delta}\\ &\propto \frac{1}{\delta^2} \end{align} ### Question (b) Show that the 95% highest posterior density region for $σ^2$ is not the same as the region obtained by squaring the endpoints of a posterior interval for σ. ### Solution 以下背理法によって示す。 c=nvとすると、2つの事後分布は以下の通り $\delta$の事後分布 \begin{align} p(\delta|data) &\propto \delta^{-1-n}exp(-c/\delta^2)\\ &= (\delta^2)^{-1/2-n/2}exp(-c/\delta^2) \end{align} $\delta^2$の事後分布 \begin{align} p(\delta|data) &\propto \delta^{-1-n}exp(-c/\delta^2)\\ &= (\delta^2)^{-1-n/2}exp(-c/\delta^2) \end{align} $(\sqrt{a},\sqrt{b})$を$p(\delta|data)$の95%区間とすると、 \begin{align} (a^2)^{-1/2-n/2}exp(-c/a) &= (b^2)^{-1/2-n/2}exp(-c/b)\\ (-1/2-n/2)\log{a}-c/a &= (-1/2-n/2)\log{b}-c/b \end{align} $(a,b)$が仮に$p(\delta^2|data)$の95%区間とすると、 \begin{align} (a^2)^{-1/2-n/2}exp(-c/a) &= (b^2)^{-1/2-n/2}exp(-c/b)\\ (-1-n/2)\log{a}-c/a &= (-1-n/2)\log{b}-c/b \end{align} 2つの式を合わせると、$1/2\log{a}=1/2\log{b}$が得られる。 したがって、$a=b$であるが、これでは区間$[a,b]$は95%区間とならず、矛盾が生じている。 ## 2.11.18 ### Question Poisson model: derive the gamma posterior distribution (2.15) for the Poisson model parameterized in terms of rate and exposure with conjugate prior distribution. ### Solution 式(2.15)(教科書p.45)ではポアソンモデルのデータ$y=y_1,...,y_i$が既知の説明変数$x=x_1,...,x_i$によって \begin{align} y_i \sim Poisson(x_i\theta) \end{align} と表されるポアソン分布に従うとき、ガンマ分布を事前分布としたときの事後分布が \begin{align} \theta|y \sim Gamma(\alpha+\sum_{i=0}^ny_i, \beta+\sum_{i=0}^nx_i) \end{align} に従うことが示されている。この導出過程を示せば良い。 $\theta \sim Gamma(\alpha, \beta)$であるから、 \begin{align} p(\theta) \propto \theta^{\alpha-1}e^{-\beta\theta} \end{align} よって、 \begin{align} p(\theta|y)&= \frac{p(y|\theta)p(\theta)}{p(y)}\\ &= \frac{Poisson(y|\theta)Gamma(\theta|\alpha,\beta)}{Poisson(y)} \\ &\propto \theta^{\sum_{i=0}^ny_i}e^{\sum_{i=0}^nx_i}\theta^{\alpha-1}e^{-\beta\theta}\\ &= \theta^{\alpha+\sum_{i=0}^ny_i}e^{-(\beta+\sum_{i=0}^nx_i)} \end{align} ## 2.11.19 ### (a) #### Question (a) Show that if y|θ is exponentially distributed with rate θ, then the gamma prior distribution is conjugate for inferences about θ given an independent and identically distributed sample of y values. #### Solution 要は尤度(y|θ ~ exp(θ))と事前分布(θ ~ Gamma(α, β))を仮定したときに事後分布がガンマ分布であることを示せば良い。$p(y|\theta)$はp56の式より \begin{align} p(y|\theta) = \theta^n e^{-n\bar{y}\theta} \end{align} θはガンマ分布であることから \begin{align} p(\theta) \propto\theta^{\alpha-1} e^{-\beta \theta} \end{align} である。このことから、 \begin{align} p(\theta|y) \propto p(y|\theta)p(\theta)\\ \propto \theta^n e^{-n\bar{y}\theta} \times \theta^{\alpha-1} e^{-\beta \theta} \\ = \theta^{n + α - 1} e^{-n\bar{y}\theta - \beta \theta}\\ = \theta^{n + α - 1} e^{-\theta(\beta + n\bar{y})} \end{align} となり、事後分布は$\theta|y \sim Gamma(α + n, β +n\bar{y})$となることが分かる。 ### (b) #### Question Show that the equivalent prior specification for the mean, φ = 1/θ, is inverse-gamma. (That is, derive the latter density function.) #### Solution (a)について、$φ = 1/θ$と定義した場合、φの事前分布が逆ガンマ分布になることを示す。 \begin{align} p(\varphi) \\ &= p(\theta)|\frac{d \varphi}{d \theta}| \\ &= \theta^{\alpha - 1} e ^{-\beta\theta} \theta^2 \\ &= \theta^{\alpha + 1} e ^{-\beta\theta} \theta^2 \\ &= \varphi^{-\alpha - 1}e^{-\frac{\beta}{\varphi}} \end{align} よって、φの事前分布は$\varphi \sim Inv-Gamma(α, β)$となることが分かる。 (解答を見てもいまいちよく理解できていないです、もし得意な方がいらっしゃれば解説をお願いしたいです) ### c #### Question The length of life of a light bulb manufactured by a certain process has an exponential distribution with unknown rate θ. Suppose the prior distribution for θ is a gamma distribution with coefficient of variation 0.5. (The coefficient of variation is defined as the standard deviation divided by the mean.) A random sample of light bulbs is to be tested and the lifetime of each obtained. If the coefficient of variation of the distribution of θ is to be reduced to 0.1, how many light bulbs need to be tested? #### Solution θの事前分布において、 \begin{align} CV(θ) \\ &= SD(θ)/E(θ) \\ &= \sqrt{Var(θ)}/E(θ)\\ &= 0.5 \end{align} である。$E(θ)$は以下のように導かれる。 \begin{align} E(θ) \\ &= \int _0 ^\infty \theta f(\theta) d\theta \\ &= \int _0 ^\infty \theta ・ \frac{\beta^\alpha}{\Gamma(\alpha)} \theta^{\alpha-1} e^{-\beta \theta} d\theta \\ &= \frac{\beta^\alpha}{\Gamma(\alpha)}\int _0 ^\infty\theta^{\alpha + 1 -1} e^{-\beta \theta} d\theta \\ &= \frac{\beta^\alpha}{\Gamma(\alpha)} \frac{\Gamma(\alpha + 1)}{\beta^{\alpha + 1}} \\ &= \frac{\alpha}{\beta} \end{align} また、(2.8)より$Var(\theta) = E((\theta)^2) - E(\theta)^2$より、$Var(\theta)$は以下の通りになる。 \begin{align} E(\theta^2) \\ &=\int_{0}^{\infty}\theta^2 \frac{\beta^\alpha}{\Gamma(\alpha)} \theta^{\alpha-1} e^{-\beta \theta}d\theta \\ &=\frac{\beta^\alpha}{\Gamma(\alpha)}\int_{0}^{\infty} \theta^{\alpha + 1} e^{-\beta \theta}d\theta \\ &= \frac{\beta^\alpha}{\Gamma(\alpha)} \frac{\Gamma(\alpha + 2)}{\beta^{\alpha +2}} \\ &= \frac{\alpha(\alpha + 1)}{\beta^2} \end{align} \begin{align} Var(\theta) \\ & = E((\theta)^2) - E(\theta)^2 \\ &= \frac{\alpha(\alpha + 1)}{\beta^2} - \frac{\alpha^2}{\beta^2}\\ &= \frac{\alpha}{\beta^2} \end{align} この2つを$CV(\theta)$に当てはめる \begin{align} CV(\theta) \\ &= \frac{\sqrt{\frac{\alpha}{\beta^2}}}{\frac{\alpha}{\beta}} \\ &= \frac{1}{\alpha^{1/2}} \\ &= 0.5 \end{align} よって$\alpha = 4$であり、aの結果より事後分布は$\theta|y \sim Gamma(α + n, β +n\bar{y})$である。よって、電球の個数をnとすると事後分布のCVは$CV(\theta|y) = \frac{1}{(\alpha + n)^{1/2}} = \frac{1}{(4 + n)^{1/2}} = 0.1$より、n = 96となる。 ### (d) #### Question In part (c), if the coefficient of variation refers to φ instead of θ, how would your answer be changed? #### Solution φの場合もcと同様のやり方でnを求めていく。逆ガンマ分布の期待値E(φ)は$E(\varphi) = \frac{\beta}{\alpha - 1}$であり、$E(\varphi^2)$は$E(\varphi^2) = \frac{\beta^2}{(\alpha - 1)(\alpha - 2)}$である。したがって  \begin{align} Var(\varphi) \\ &= E(\varphi^2) - E(\varphi)^2 \\ &= \frac{\beta^2}{(\alpha - 1)(\alpha - 2)} - \frac{\beta^2}{(\alpha - 1)^2} \\ &= \frac{\beta^2}{(\alpha - 1)^2(\alpha - 2)} \end{align} よって、CVは$\frac{\sqrt(\frac{\beta^2}{(\alpha - 1)^2(\alpha - 2)})}{\frac{\beta}{\alpha - 1}} = (\alpha -2)^{-1/2} = 0.5$になるから、$\alpha = 6$。 事後分布のCVは、$(\alpha + n + 2 - 2)^{-1/2} = 0.1$であることから、n = 94となる。 ## 2.11.20 ### (a) #### Question Suppose y|θ is exponentially distributed with rate θ, and the marginal (prior) distri- bution of θ is Gamma(α,β). Suppose we observe that y ≥ 100, but do not observe the exact value of y. What is the posterior distribution, p(θ|y≥100), as a function of α and β? Write down the posterior mean and variance of θ. #### Solution $y|\theta \sim Exp(\theta)$と$\theta \sim Gamma(\alpha,\beta)$は前問と同様である。 $p(y|\theta) = \theta e^{-y\theta}$,$p(\theta) \propto\theta^{\alpha-1} e^{-\beta \theta}$より、 \begin{align} p(\theta|y \geq 100) \\ & \propto\theta^{\alpha-1} e^{-\beta \theta}\int_{100}^{\infty} \theta e^{-y\theta}dy \\ &= \theta^{\alpha-1} e^{-\beta \theta}\left[\frac{1}{-1}e^{-y\theta}\right]^\infty_{100} \\ &= \theta^{\alpha-1}e^{-\beta \theta}e^{-100\theta} \\ &= \theta^{\alpha-1}e^{-(\beta + 100) \theta} \end{align} よって、事後分布は$\theta|y \sim Gamma(α, β + 100)$となる。前問より、平均と分散はそれぞれ$\frac{\alpha}{\beta + 100}$, $\frac{\alpha}{(\beta + 100)^2}$になる。 ### (b) #### Question In the above problem, suppose that we are now told that y is exactly 100. Now what are the posterior mean and variance of θ? #### Solution aについてy = 100の場合を考える。 \begin{align} p(\theta|y = 100) \\ &\propto p(y = 100|\theta)p(\theta)\\ &= \theta e^{-100\theta} \theta^{\alpha -1}e^{-\beta\theta} \\ &= \theta^\alpha e^{-\theta(\beta + 100)} \end{align} よって、y = 100の時事後分布は$\theta|y \sim Gamma(α + 1, β + 100)$となる。前問より、平均と分散はそれぞれ$\frac{\alpha + 1}{\beta + 100}$, $\frac{\alpha + 1}{(\beta + 100)^2}$になる。 ### (c) #### Question Explain why the posterior variance of θ is higher in part (b) even though more in- formation has been observed. Why does this not contradict identity (2.8) on page 32? #### Solution 確かに、aの事後分散$\frac{\alpha}{(\beta + 100)^2}$はbの事後分散$\frac{\alpha + 1}{(\beta + 100)^2}$よりも小さい。それは、(2.8)から示された$Var(\theta|y \geqq 100)\geqq E(Var(\theta|y)| y \geqq 100)$と矛盾しているように見える。しかし、平均的には上記式の右辺の方が小さくなると言っているだけであり、ある1つの値についても同様のことが言えるわけではない。

    Import from clipboard

    Paste your markdown or webpage here...

    Advanced permission required

    Your current role can only read. Ask the system administrator to acquire write and comment permission.

    This team is disabled

    Sorry, this team is disabled. You can't edit this note.

    This note is locked

    Sorry, only owner can edit this note.

    Reach the limit

    Sorry, you've reached the max length this note can be.
    Please reduce the content or divide it to more notes, thank you!

    Import from Gist

    Import from Snippet

    or

    Export to Snippet

    Are you sure?

    Do you really want to delete this note?
    All users will lose their connection.

    Create a note from template

    Create a note from template

    Oops...
    This template has been removed or transferred.
    Upgrade
    All
    • All
    • Team
    No template.

    Create a template

    Upgrade

    Delete template

    Do you really want to delete this template?
    Turn this template into a regular note and keep its content, versions, and comments.

    This page need refresh

    You have an incompatible client version.
    Refresh to update.
    New version available!
    See releases notes here
    Refresh to enjoy new features.
    Your user state has changed.
    Refresh to load new user state.

    Sign in

    Forgot password

    or

    By clicking below, you agree to our terms of service.

    Sign in via Facebook Sign in via Twitter Sign in via GitHub Sign in via Dropbox Sign in with Wallet
    Wallet ( )
    Connect another wallet

    New to HackMD? Sign up

    Help

    • English
    • 中文
    • Français
    • Deutsch
    • 日本語
    • Español
    • Català
    • Ελληνικά
    • Português
    • italiano
    • Türkçe
    • Русский
    • Nederlands
    • hrvatski jezik
    • język polski
    • Українська
    • हिन्दी
    • svenska
    • Esperanto
    • dansk

    Documents

    Help & Tutorial

    How to use Book mode

    Slide Example

    API Docs

    Edit in VSCode

    Install browser extension

    Contacts

    Feedback

    Discord

    Send us email

    Resources

    Releases

    Pricing

    Blog

    Policy

    Terms

    Privacy

    Cheatsheet

    Syntax Example Reference
    # Header Header 基本排版
    - Unordered List
    • Unordered List
    1. Ordered List
    1. Ordered List
    - [ ] Todo List
    • Todo List
    > Blockquote
    Blockquote
    **Bold font** Bold font
    *Italics font* Italics font
    ~~Strikethrough~~ Strikethrough
    19^th^ 19th
    H~2~O H2O
    ++Inserted text++ Inserted text
    ==Marked text== Marked text
    [link text](https:// "title") Link
    ![image alt](https:// "title") Image
    `Code` Code 在筆記中貼入程式碼
    ```javascript
    var i = 0;
    ```
    var i = 0;
    :smile: :smile: Emoji list
    {%youtube youtube_id %} Externals
    $L^aT_eX$ LaTeX
    :::info
    This is a alert area.
    :::

    This is a alert area.

    Versions and GitHub Sync
    Get Full History Access

    • Edit version name
    • Delete

    revision author avatar     named on  

    More Less

    Note content is identical to the latest version.
    Compare
      Choose a version
      No search result
      Version not found
    Sign in to link this note to GitHub
    Learn more
    This note is not linked with GitHub
     

    Feedback

    Submission failed, please try again

    Thanks for your support.

    On a scale of 0-10, how likely is it that you would recommend HackMD to your friends, family or business associates?

    Please give us some advice and help us improve HackMD.

     

    Thanks for your feedback

    Remove version name

    Do you want to remove this version name and description?

    Transfer ownership

    Transfer to
      Warning: is a public team. If you transfer note to this team, everyone on the web can find and read this note.

        Link with GitHub

        Please authorize HackMD on GitHub
        • Please sign in to GitHub and install the HackMD app on your GitHub repo.
        • HackMD links with GitHub through a GitHub App. You can choose which repo to install our App.
        Learn more  Sign in to GitHub

        Push the note to GitHub Push to GitHub Pull a file from GitHub

          Authorize again
         

        Choose which file to push to

        Select repo
        Refresh Authorize more repos
        Select branch
        Select file
        Select branch
        Choose version(s) to push
        • Save a new version and push
        • Choose from existing versions
        Include title and tags
        Available push count

        Pull from GitHub

         
        File from GitHub
        File from HackMD

        GitHub Link Settings

        File linked

        Linked by
        File path
        Last synced branch
        Available push count

        Danger Zone

        Unlink
        You will no longer receive notification when GitHub file changes after unlink.

        Syncing

        Push failed

        Push successfully