36-711 Statistical Computing Homework 4 Due Nov XX, 2002

1. If the complete data X follows the exponential family, then the likelihood for X is of

the form

f(X|θ) = A(X)e^{B(θ)+K(X)^{T} γ(θ)}

Suppose Y = t(X) is the observed data, and define

μ(θ) = E[K(X)|θ]  complete-data expected suff. stats
μ_{k}(Y) = E[K(X)|Y, θ^{k}] missing-data expected suff. stats

(a) Show that the E-M algorithm for estimating θ in this exponential family may be expressed as:

Setup: Compute μ(θ) and an initial guess θ_{0}

For each step k:

E step: Compute μ_{k}(Y)
M step: θ^{k+1} = argmax_{θ}B(θ) + μ_{k}(Y)^{T} γ(θ)

(b) Use the fact that B(θ) is the log-normalizing constant for the model f(X|θ), to

show that

∇B(θ) = −μ^{T} (θ)∇γ(θ)

Hence, show that setting the gradient equal to zero to find θ^{k+1} in the M step leads to the following restatement of the E-M algorithm:

Setup: Compute μ(θ) and an initial guess θ_{0}

For each step k:

E step: Compute μ_{k}(Y)
M step: θ^{k+1} = root_{θ} { μ^{T} (θ)∇γ(θ) = μ_{k}^{T} (Y)∇γ(θ) }

If θ are the natural parameters for the model, then γ(θ) = θ. In that case, the M step is just setting θ^{k+1} = root_{θ} { μ(θ) = μ_{k}(Y) }. How does this compare with maximum likelihood in exponential families with no missing data?

2. In a Bayesian model for complete data X there will be a likelihood f(X|θ) and a prior distribution π(θ), and then the complete-data posterior distribution for θ will be of the form

π(θ|X) ∝ f(X|θ)π(θ)

(a) Show that if Y = t(X) is the observed data, then the observed-data posterior for

θ is of the form

π(θ|Y) ∝ g(Y |θ)π(θ)

where g(Y |θ) is the complete-data likelihood.

(b) Define Q_{post}(θ, θ^{k}) = E[π(X|θ)|Y, θ^{k}], and show (directly, or by analogy with or extension of the likelihood case considered in class) that

log π(θ|Y) − Q_{post}(θ, θ^{k})

achieves its minimum at θ = θ^{k}. Hence the E-M algorithm

E step: Compute Q_{post}(θ, θ^{k})
M step: θ^{k+1} = argmax_{θ}Q_{post}(θ, θ^{k})

finds a posterior mode (or saddle point) θ̂_{post} of the observed-data posterior π(θ|Y).

(c) The quantities −∇^{2} log π(θ|Y) and I_{post}(θ) = E[−∇^{2} log π(θ|Y)|Y] function as “observed” and “expected” posterior information functions and can be inverted to produce posterior standard errors for θ̂_{post} [the proof of this involves the same sort of taylor expansion argument that leads to asymptotic normality of MLE’s]. Show, by analogy with or extension of the likelihood case discussed in class, that the following Missing Information Principle holds:

−∇^{2} log π(θ|Y) = E[−∇^{2} log π(θ|X)|Y, θ] − Var[∇ log π(θ|X)|Y, θ]

3. Let t_{1}, …, t_{n} be iid with f(t) = ρe^{−ρt}, the exponential distribution with parameter ρ. Suppose that we observe (y_{i}, δ_{i}) where y_{i} = min(t_{i}, c_{i}) and δ_{i} = 1 if t_{i} < c_{i} (t_{i} is uncensored), and δ_{i} = 0 otherwise (t_{i} is censored), for i = 1, …, n. We want to use E-M to estimate ρ via maximum likelihood.

(a) Show that the E-step computes

Q(ρ, ρ^{k}) = n log ρ − ρ ∑_{i−1}^{n} E[t_{i}|y_{i}, δ_{i}, ρ^{k}]

(b) For an uncensored observation, t_{i} = y_{i}. What is E[t_{i}|y_{i}, δ_{i}, ρ^{k}] for a censored observation?

(c) What is the M-step for this problem?

(d) Implement the E-M algorithm using the data given below (survivial time after transplant in the Stanford Heart Transplant study). How did you assess convergence? How many iterations were required? Try several different or randomly selected starting values. Is there evidence that the observed likelihood is multimodal?

(e) Compute a standard error for the best ρ̂ you can come up with, using the missing information principle.

(f) The full data set has 184 observations. Compare your results for the first 45 observations (about 1/4 of the data) with your results for the full data set.

Stanford Heart Transplant Data

linux% cp ~/brian/F711/stan.R .
linux% Splus   # or R...
> source("stan.R")
> stan[1:5,]
  time status age
1   15      1  54
2    3      1  40
3   46      1  42
4  623      1  51
5  126      1  48
> dim(stan)
[1] 184   3

“time” is y_{i}, “status” (dead=1, alive=0) is δ_{i}.

4. Saying

> library(MASS)     # may skip this step in R
> data(faithful)    # skip this step in Splus
> attach(faithful)

in R or Splus defines a vector waiting of n = 272 waiting times between eruptions of the geyser Old Faithful at Yellowstone Park. [A related data set, geyser is discussed by Venables and Ripley (MASS 3rd ed, pp 252ff; 4th ed. pp. 127ff)].

(a) Make a histogram of the data in waiting. You will see that a plausible model for the waiting times y_{1}, …, y_{n} is a mixture of two normals:

\begin{aligned}y_{i} | μ_{1}, μ_{2}, σ_{1}, σ_{2}, q & ∼^{iid} & \{q \frac{1}{\sqrt{2π}σ_{1}} exp[−(y_{i} − μ_{1})^{2}/2σ_{1}^{2}] \\ & & + (1 − q)\frac{1}{\sqrt{2π}σ_{2}} exp[−(y_{i} − μ_{2})^{2}/2σ_{2}^{2}]\} \end{aligned}  (∗)

We wish to estimate the parameters μ_{1}, μ_{2}, σ_{1}, σ_{2} and q. Explicitly write the likelihood g(\underline{y}|μ_{1}, μ_{2}, σ_{1}, σ_{2}, q) for the observed data \underline{y} = (y_{1}, …, y_{n}).

(b) The likelihood is greatly simplified by using data augmentation. Consider the latent variable (a.k.a. missing data, or augmentation data)

z_{i} = \{ \begin{matrix}1 & if y_{i} belongs to the left hump \\ 0 & if y_{i} belongs to the right hump \end{matrix} , i = 1, …, n .

Show that if we assume

z_{i} | q ∼^{iid} Bin(1, q)
y_{i} | μ_{1}, μ_{2}, σ_{1}, σ_{2}, z_{i} ∼^{indep} \{ \begin{matrix}N(μ_{1}, σ_{1}^{2}), & if z_{i} = 1 \\ N(μ_{2}, σ_{2}^{2}), & if z_{i} = 0 \end{matrix}

then the distribution of the observed data \underline{y} = (y_{1}, …, y_{n}) is as in (∗). Show, therefore, that we may take the complete-data likelihood for the data vectors \underline{y} and \underline{z} to be

f(\underline{y}, \underline{z}|μ_{1}, μ_{2}, σ_{1}, σ_{2}, q) = ∏_{i=1}^{n} \{q \frac{1}{\sqrt{2π}σ_{1}} exp[−(y_{i} − μ_{1})^{2}/2σ_{1}^{2}]\}^{z_{i}}

× \{(1 − q)\frac{1}{\sqrt{2π}σ_{2}} exp[−(y_{i} − μ_{2})^{2}/2σ_{2}^{2}]\}^{1−z_{i}}

(c) Derive the function E-M objective function Q(θ, θ^{k}) for this problem, where θ = (μ_{1}, μ_{2}, σ_{1}, σ_{2}, q). Are there any simplifications to the E-step? The M-step? Notes:

• z_{i} | y_{i}, μ_{1}, μ_{2}, σ_{1}, σ_{2}, q are still independent Bernoulli’s but now with different success probabilities; what are they?

• A modification of E-step that might be useful is to round the E[z_{i}|y_{i}, μ_{1}, μ_{2}, σ_{1}, σ_{2}, q] to 0 or 1 rather than use the expected value itself. How does this simplify the M-step?

(d) Fit the model to the Old Faithful data, using E-M or some variation. Describe your findings, focusing on the following questions:

• What would be a good set of starting values?

• How do you assess convergence?

• Try from several different sets of staring values. Do you get different maxima? Which maximum or maxima are “best”? Why?

(e) Obtain standard errors for the parameters two ways: using infomat.s, and using the missing information principle.