Poisson population comparisons : Let’s reconsider the number of children data of Exercise 4.8. We’ll assume Poisson sampling models for the two groups as before, but now we’ll parameterize $θ_A$ and $θ_B$ as $θ_A = θ$, $θ_B = θ × γ.$ In this parameterization, $γ$ represents the relative rate $θ_B/θ_A$. Let $θ ∼ gamma(a_θ, b_θ)$ and let $γ ∼ gamma(a_γ, b_γ)$.
Group A
Average number of children in their 30s with a bachelor’s degree
$\sim Poission(\theta_A)=Poisssion(\theta)$
Group B
Average number of children in their 30s without a bachelor’s degree
$\sim Poission(\theta_B)=Poisssion(\theta\cdot \gamma)$
$\theta \sim gamma(a_\theta, b_\theta)$
$\gamma \sim gamma(a_\gamma, b_\gamma)$
$\theta, \gamma$ is independent
a) Are $θ_A$ and $θ_B$ independent or dependent under this prior distribution? In what situations is such a joint prior distribution justified?
i) Are $\theta_A$ and $\theta_B$ independent?
If $\theta_A$ and $\theta_B$ are independent, covariance = 0
$Cov(\theta_A , \theta_B) = E[\theta_A \theta_B]-E[\theta_A]E[\theta_B]$
$= E[\theta \cdot \theta \gamma]-E[\theta]E[\theta\cdot\gamma]$
$= E[\theta^2]E[\gamma]-E[\theta]E[\theta]E[\gamma] \;\;\;\;\because \theta \;\; and \;\; \gamma \;\; are \;\; independent$
$= E[\theta^2]E[\gamma]-E[\theta]^2E[\gamma]$
$=(E[\theta^2]-E[\theta]^2)E[\gamma]$
$=Var[\theta] \cdot E[\gamma]$
$\therefore Cov(\theta_A,\theta_B) = Var(\theta) \cdot E[\gamma]=a_\theta b^2_{\theta}\cdot a_\gamma b_\gamma \neq0$
$\rightarrow \theta_A$ and $\theta_B$ are not independent
ii)
b) Obtain the form of the full conditional distribution of $θ$ given $y_A, y_B$ and $γ$.
$p(\theta|y_A,y_B,\gamma) \propto \;p(y_A,y_A|\theta,\gamma)\times p(\theta,\gamma)$
$= p(y_A |\theta)p(y_B|\theta,\gamma)\times p(\theta) \times p( \gamma) \; \because y_A\;and\; \gamma \; are \; independent$
$\propto (\prod_{i=1}^{n_A} \theta^{y_{A_i}}e^{-\theta})\;(\prod_{i=1}^{n_B} (\theta\gamma)^{y_{B_i}}e^{-\theta\gamma}) \times (\theta^{a_{\theta}-1}e^{-b_{\theta}\theta})\times(\gamma^{a_r-1}e^{-b_\gamma \gamma})$
$=(\theta^{a_\theta-1}e^{-b_\theta \theta})\times (\gamma^{a_\gamma-1}e^{b_\gamma \gamma})\times(\theta^{\sum^{n_A}_{i=1}y_{A_i}e^{-n_A \theta}})\times (\gamma^{\sum^{n_B}_{i=1}y_{B_i}e^{-n_B \gamma}})$
$=(\theta^{a_\theta-1}e^{-b_\theta \theta})\times (\gamma^{a_\gamma-1}e^{b_\gamma \gamma})\times(\theta^{n_a \bar{y}_A}e^{-n_A\theta})\times ((\gamma \theta)^{n_B \bar{y}_B}e^{-n_B \gamma \theta})$
$\propto(\theta^{a_\theta-1}e^{-b_\theta \theta})\times(\theta^{n_a \bar{y}_A}e^{-n_A\theta})\times ((\gamma \theta)^{n_B \bar{y}_B}e^{-n_B \gamma \theta})$
$\propto \theta^{(a_{\theta}+n_A\bar{y}_{A}+n_B\bar{y}_{B})-1}\;e^{-(n_A+n_B\gamma+b_{\theta})\theta}$
$\therefore p(\theta|y_A,y_B,\gamma) \propto dgamma(a_\theta+n_A\bar{y}+n_B\bar{y}B, b\theta+n_A+n_B\gamma)$
c) Obtain the form of the full conditional distribution of $γ$ given $y_A, y_B$ and $θ$.
$p(\gamma\;|\;y_A,y_B,\theta) \propto p(y_A,y_A|\theta,\gamma)\times p(\theta,\gamma)$
$\propto p(y_A |\theta)p(y_B|\theta,\gamma)\times p(\theta) \times p( \gamma) \; \because y_A\;and\; \gamma \; are \; independent$
$\propto (\prod_{i=1}^{n_A} \theta^{y_{A_i}}e^{-\theta})\;(\prod_{i=1}^{n_B} (\theta\gamma)^{y_{B_i}}e^{-\theta\gamma}) \times (\theta^{a_{\theta}-1}e^{-b_{\theta}\theta})\times(\gamma^{a_r-1}e^{-b_\gamma \gamma})$
$\propto (\gamma^{a_\gamma -1}e^{-b_\gamma \gamma})\times (\gamma^{n_B\bar{y}_B}e^{-n_B\gamma_\theta})$
$\propto \gamma^{(a_{\gamma}+n_B\bar{y}_{B})-1}\;e^{-(n_B\theta+b_{\gamma})\gamma}$
$\propto dgamma(a_\gamma+n_B\bar{y}_B,n_B\theta+b_\gamma)$
d) Set $a_θ = 2$ and $b_θ = 1$. Let $a_γ = b_γ ∈ \{8, 16, 32, 64, 128\}$. For each of these five values, run a Gibbs sampler of at least 5,000 iterations and obtain $E[θ_B −θ_A|y_A, y_B]$. Describe the effects of the prior distribution for $γ$ on the results.
y_a <- scan('<https://www2.stat.duke.edu/courses/Fall09/sta290/datasets/Hoffdata/menchild30bach.dat>')
y_b <- scan('<https://www2.stat.duke.edu/courses/Fall09/sta290/datasets/Hoffdata/menchild30nobach.dat>')
len_a <- length(y_a)
len_b <- length(y_b)
bar_y_a <- mean(y_a)
bar_y_b <- mean(y_b)
theta_a <- 2
theta_b <- 1
sample <- 5000
gamma_ab <- c(8, 16, 32, 64, 128)
theta <- c() ; gamma <- c()
initial_theta <- bar_y_a
initial_gamma <- bar_y_a / bar_y_b
gibbs_sampler <- function(gamma_ab_value){
for (s in 1:sample) {
gamma_a <- gamma_ab_value; gamma_b <- gamma_ab_value
initial_theta = rgamma(
1, theta_a + len_a * bar_y_a + len_b * bar_y_b, theta_b + len_a + len_b * initial_gamma
)
initial_gamma = rgamma(
1, gamma_a + len_b * bar_y_b, gamma_b + len_b * initial_theta
)
theta[s] <- initial_theta
gamma[s] <- initial_gamma
theta_A <- theta
theta_B <- theta * gamma
}
mean_ab <- mean(theta_B - theta_A)
return(mean_ab)
}
for (i in gamma_ab) {
gibbs <- gibbs_sampler(i)
means <- print(gibbs)
means
}
means <- c(gibbs_sampler(8), gibbs_sampler(16), gibbs_sampler(32), gibbs_sampler(64), gibbs_sampler(128))
plot(gamma_ab, means, type='o')
axis(1, at= gamma_ab)
Mixture model: The file glucose.dat contains the plasma glucose concentration of 532 females from a study on diabetes (see Exercise 7.6).
a) Make a histogram or kernel density estimate of the data. Describe how this empirical distribution deviates from the shape of a normal distribution.
It’s possible to figure out that data skewed right.
b) Consider the following mixture model for these data: For each study participant there is an unobserved group membership variable $X_i$ which is equal to 1 or 2 with probability $p$ and $1 − p$.
If $X_i = 1$ then $Y_i ∼ normal(θ_1, σ^2_1)$, and if $X_i = 2$ then $Y_i ∼ normal(θ_2, σ^2_2)$.
Let $p ∼ beta(a, b)$, $θ_j ∼ normal(µ_0, τ^2_0)$ and $1/σ_j ∼ gamma(ν_0/2, ν_0σ^2_0/2)$ for both $j = 1$ and $j = 2$.
Obtain the full conditional distributions of $(X_1, . . . , X_n), p, θ_1, θ_2, σ^2_1$ and $σ^2_2$.
The likelihood is
$p(y|x,p,\theta_1,\theta_2,\sigma_1^2, \sigma_2^2)= \prod_{i=1}^n p(y_i|x_i,p,\theta_1,\theta_2,\sigma_1^2, \sigma_2^2)$
$=\prod_{i=1}^n [dnorm(y_i, \theta_1,\sigma_1^2)]^{x_i} \times [dnorm(y_i, \theta_2,\sigma_2^2)]^{1-x_i}$
Full conditional distribution
$p(\theta_1 \;|\;y, \theta_2) \;\propto\;p(y \;|\;\theta_1,\theta_2)\times p(\theta_1)$
i) In case of a single $X_i =1$
$P(X_i = 1|y_i, p, \theta_1, \theta_2, \sigma^2_1, \sigma_2^2)=\frac{P(X_i=1|p, \theta_1, \theta_2, \sigma_1^2, \sigma_2^2)\times P(y_i|X_i=1,p, \theta_1, \theta_2, \sigma_1^2, \sigma_2^2) }{P(y_i|p, \theta_1, \theta_2, \sigma_1^2, \sigma_2^2)}$
$= \frac{P(X_1=1|p)\times p(y_i|X_1=1, \theta_1, \sigma^2_1)}{P(X_i=1|p)\times p(y_i|X_i=1, \theta_1, \sigma^2_1)+P(X_i=0|p)\times p(y_i|X_i=0, \theta_2, \sigma^2_2)}$
$= \frac{p \times dnorm(y_i, \theta_1, \sigma^2_2)}{p\times dnorm(y_i, \theta_1, \sigma^2_1)+(1-p)\times dnorm(y_i, \theta_2, \sigma^2_2)}$
$P(X_i=2|y_i, p, \theta_1, \theta_2, \sigma^2_1, \sigma^2_2)=$ ${p(X_i=2\;|\;p,\theta_1,\theta_2,\sigma_1^2, \sigma_2^2)\;\;\times\;\; p(y_i\;|\;X_i=2,p,\theta_1,\theta_2,\sigma_1^2, \sigma_2^2) \over p(y_i\;|\;p,\theta_1,\theta_2,\sigma_1^2, \sigma_2^2)}$
$={p(X_i=2\;|\;p)\;\;\times\;\; p(y_i\;|\;X_i=2,\theta_1,\sigma_1^2) \over p(X_i=2\;|\;p)\times p(y_i\;|\;X_i=2,\theta_1,\sigma_1^2) \;+\; p(X_i=1\;|\;p)\times p(y_i\;|\;X_i=2, \theta_2,\sigma_2^2)} \\\\$
$=
{p\times dnorm(y_i, \theta_1,\sigma_1^2) \over p\times dnorm(y_i, \theta_1,\sigma_1^2) \;+\; (1-p)\times dnorm(y_i, \theta_2,\sigma_2^2)}$
$\because$ it’s similar to $P(X_i=0)$
$x_i \sim$ $p(x_i \;|\; y_i,p,\theta_1,\theta_2,\sigma_1^2, \sigma_2^2) \;\propto\;Bern({p \times dnorm(y_i,\theta_1,\sigma_1^2) \over p\times dnorm(y_i,\theta_1,\sigma_1^2)\;+\;(1-p)\times dnorm(y_i,\theta_2,\sigma_2^2)})$
Let if $n_1=\sum x_i$, and $n_2=n-n_1$,
$p(p\;|\;x,y,\theta_1,\theta_2,\sigma_1^2, \sigma_2^2)\propto p(p)\times p(x,y,\theta_1,\theta_2,\sigma_1^2, \sigma_2^2|p)$
$= p(p)\times p(x|p)p(y|x,\theta_1,\theta_2,\sigma_1^2, \sigma_2^2)p(\theta_1,\theta_2,\sigma_1^2, \sigma_2^2) \\\\$
$=
p(p)\times p(x|p)$
$= dbeta(p,a,b) \times dbinom(\sum x_i,n,p) \\\\$
$= p^{a-1}(1-p)^{b-1} \times p^{\sum x_i}(1-p)^{n-\sum x_i} \\\\$
$= p^{(a+\sum x_i)-1}(1-p)^{(b+n-\sum x_i)-1} \\\\$
$= dbeta(p,\;\;a+\sum x_i,\;\;b+n-\sum x_i)$
$\therefore p(p|x,y,\theta_1,\theta_2,\sigma_1^2, \sigma_2^2) \;\;\propto\;\;dbeta(p,\;\;a+\sum x_i,\;\;b+n-\sum x_i)$
Let $y_i =$ {$y_i \in y, x_i = 1$} and $y_2=${$y_i \in y, x_i = 0$}
where, $v_{n,1}=v_0, n_1$ and $\sigma^2_{n,1}(\theta_1)=\frac{1}{v_{n,1}}[v_0 \sigma^2_0 + n_1 s^2_{n,1}(\theta_1)]$,
$p(\theta_1|x,y,p,\theta_2,\sigma_1^2, \sigma_2^2) \propto p(\theta_1|x,p,\theta_2,\sigma_1^2, \sigma_2^2)\times p(y|x,p,\theta_2,\sigma_1^2, \sigma_2^2)$
$= p(\theta_1)\times \prod_{i=1}^n dnorm(y_i,\theta_1,\sigma_1^2)$
$= dnorm(\theta_1, \mu_0, \tau_0^2) \times \prod_{i=1}^n dnorm(y_i,\theta_1,\sigma_1^2)$
$= exp(-{1\over \tau_0^2}(\theta_1 - \mu_0)^2) \times \prod_{y \in y_i} exp(-{1\over2\sigma_1^2\sum(y_i-\theta_1)^2})$
$\sim Norm(\;1/({1\over\tau_0^2}+{\sum x_i \over \sigma_1^2})\;,\;({1\over\tau_0^2}\mu_0 + {\sum x_i\over\sigma_1^2}\bar{y})/({1\over\tau_0^2}+{\sum x_i \over \sigma_1^2}))$
$\propto Norm(\tau_n^2, \mu_{n,1})$
When we switch from 1 to 2, it’s possible to get $\theta_2 | \cdots Norm(\mu_{n,2}, \tau^2_{n,2)}$
$p(\sigma_1^2|x,y,p, \theta_1, \theta_2, \sigma^2_2) \propto p(\sigma^2_1|x, p, \theta_1, \theta_2, \sigma^2_2)\times p(y|x, p, \theta_1, \theta_2, \sigma^2_1, \sigma^2_2)$
$\propto p(\sigma^2_1) \times \prod_{i=1}^n(dnorm(y_i, \theta_1, \sigma^2_1)^{x_i} dnorm(y_i, \theta_2, \sigma^2_2)^{1-x_i})$
$\propto IG(\sigma_1^2,\nu_0,\nu_0\sigma_0^2/2) \times \prod_{i=1}^n dnorm(y_i,\theta_1,\sigma_1^2) \\\\$
$= exp((\sigma_1^2)^{-(\nu_0/2)-1} \;exp(-{1\over\sigma_1^2}\nu_0\sigma_0^2/2))$
$\times (\sigma_1^2)^{-n/2} exp(-{1\over2\sigma_1^2}\sum (y_i-\theta_1)^2)$
$\propto IG(v_{n,1}/2, \sigma^2_{n, 1}(\theta_1)v_{n,1}/2$