\documentstyle[12pt,fullpage]{article}

\let\grad\nabla
\let\del\partial

\def\Var{\mbox{Var}}
\def\argmax{\mbox{argmax}}

\def\beq{$$}
\def\eeq{$$}
\def\bea{\begin{eqnarray*}}
\def\eea{\end{eqnarray*}}
\let\ul\underline

\begin{document}

\parskip=0.6em
\parindent=0pt

\begin{center}
\large\bf
36-711 Statistical Computing \\
Homework 4 \\
Due Nov XX, 2002
\end{center}

\begin{enumerate}

\item If the complete data $X$ follows the exponential family, then
  the likelihood for $X$ is of the form
\[
        f(X|\theta) = A(X) e^{B(\theta) + K(X)^T\gamma(\theta)}
\]
Suppose $Y=t(X)$ is the observed data, and define
\begin{eqnarray*}
        \mu(\theta) & = & E[K(X) | \theta]  ~~~~~~~~~ \mbox{complete-data
        expected suff. stats} \\
        \mu_k(Y)    & = & E[K(X) | Y, \theta^k] ~~~~ \mbox{missing-data
        expected suff. stats}  
\end{eqnarray*}
\begin{enumerate}
\item Show that the E-M algorithm for estimating $\theta$ in this
        exponential family may be expressed as:
\begin{eqnarray*}
\mbox{Setup:}  & & \mbox{Compute $\mu(\theta)$ and an initial guess
  $\theta_0$}\\
\mbox{For each step $k$:} & & \\
\mbox{E step:} & & \mbox{Compute } \mu_k(Y)\\
\mbox{M step:} & & \theta^{k+1} = \argmax_{\theta} B(\theta) +
  \mu_k(Y)^T\gamma(\theta) 
\end{eqnarray*}

\item Use the fact that $B(\theta)$ is the log-normalizing constant
  for the model $f(X|\theta)$, to show that
\[
        \grad B(\theta) = - \mu^T(\theta)\grad\gamma(\theta)
\]
Hence, show that setting the gradient equal to zero to find
$\theta^{k+1}$ in the M step leads to the following restatement of the
E-M algorithm:
\begin{eqnarray*}
\mbox{Setup:}  & & \mbox{Compute $\mu(\theta)$ and an initial guess
  $\theta_0$}\\
\mbox{For each step $k$:} & & \\
\mbox{E step:} & & \mbox{Compute } \mu_k(Y)\\
\mbox{M step:} & & \theta^{k+1} = \mbox{root}_\theta\left\{~
         \mu^T(\theta)\grad\gamma(\theta) ~ = ~  
            \mu_k^T(Y)\grad\gamma(\theta)
          ~\right\}
\end{eqnarray*}
If $\theta$ are the natural parameters for the model, then
$\gamma(\theta) = \theta$.  In that case, the M step is just setting
$\theta^{k+1} = \mbox{root}_\theta\left\{~ \mu(\theta) ~=~
\mu_k(Y)~\right\}$.  How does this compare with maximum likelihood in
exponential families with no missing data?

\end{enumerate}

\item In a Bayesian model for complete data $X$ there will be a
  likelihood $f(X|\theta)$ and a prior distribution $\pi(\theta)$, and
  then the complete-data posterior distribution for $\theta$ will be
  of the form 
\[
        \pi(\theta|X) \propto f(X|\theta) \pi(\theta)
\]

\begin{enumerate}

\item Show that if $Y=t(X)$ is the observed data, then the
  observed-data posterior for $\theta$ is of the form
\[
        \pi(\theta|Y) \propto g(Y|\theta) \pi(\theta)
\]
where $g(Y|\theta)$ is the complete-data likelihood.

\item Define $Q_{post}(\theta,\theta^k) =
  E[\pi(X|\theta)|Y,\theta^k]$, and show (directly, or by analogy with
  or extension of the likelihood case considered in class) that
\[
        \log\pi(\theta|Y) - Q_{post}(\theta,\theta^k)
\]
achieves its minimum at $\theta=\theta^k$.  Hence the E-M algorithm
\begin{eqnarray*}
\mbox{E step:} & & \mbox{Compute } Q_{post}(\theta,\theta^k)\\
\mbox{M step:} & & \theta^{k+1} = \argmax_{\theta} Q_{post}(\theta,\theta^k)
\end{eqnarray*}
finds a posterior mode (or saddle point) $\hat\theta_{post}$ of
the observed-data posterior $\pi(\theta|Y)$. 

\item The quantities $-\grad^2 \log\pi(\theta|Y)$ and $I_{post}(\theta) =
  E[-\grad^2\log\pi(\theta|Y)|Y]$ function as ``observed'' and
  ``expected'' posterior information functions and can be inverted to
  produce posterior standard errors for $\hat\theta_{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:
\[
        -\grad^2\log\pi(\theta|Y) = 
        E[-\grad^2\log\pi(\theta|X)|Y,\theta] -
        \Var[\grad\log\pi(\theta|X)|Y,\theta]
\]

\end{enumerate}


\item Let $t_1, \ldots, t_n$ be iid with  $f(t) =
  \rho e^{-\rho t}$, the exponential distribution with parameter
  $\rho$.  Suppose that we observe $(y_i,\delta_i)$ where
  $y_i=min(t_i,c_i)$ and $\delta_i=1$ if $t_i<c_i$ ($t_i$ is {\em
    uncensored}), and $\delta_i=0$ otherwise ($t_i$ is {\em
    censored}), for $i=1,\ldots, n$.  We want to use E-M to estimate
  $\rho$ via maximum likelihood.
\begin{enumerate}

\item Show that the E-step computes
\[
        Q(\rho,\rho^k) = n\log \rho - \rho \sum_{i-1}^n 
        E[t_i|y_i,\delta_i,\rho^k]
\]

\item For an uncensored observation, $t_i=y_i$.  What is
        $E[t_i|y_i,\delta_i,\rho^k]$ for a censored observation?

\item What is the M-step for this problem?

\item 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?

\item Compute a standard error for the best $\hat\rho$  you can come
  up with, using the missing information principle.

\item 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.

\end{enumerate}

\begin{center}
\begin{tabular}{|p{4in}|}
\hline
Stanford Heart Transplant Data
\begin{verbatim}
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
\end{verbatim}
``time'' is $y_i$, ``status'' (dead=1, alive=0) is $\delta_i$.  
\\ \hline
\end{tabular}
\end{center}




\item Saying
\begin{verbatim}
> library(MASS)     # may skip this step in R
> data(faithful)    # skip this step in Splus
> attach(faithful)
\end{verbatim}
  in R or Splus defines a vector {\tt waiting} of $n=272$ waiting times
  between eruptions of the geyser Old Faithful at Yellowstone Park.
  [A related data set, {\tt geyser} is discussed by Venables and
  Ripley (MASS 3rd ed, pp 252ff; 4th ed.  pp. 127ff)].

\begin{enumerate}

\item Make a histogram of the data in {\tt waiting}.
  You will see that a plausible model for the waiting times $y_1,
\ldots, y_{n}$ is a mixture of two normals:
\beq
        \begin{array}{rcl}
         y_i ~|~\mu_1, \mu_2, \sigma_1, \sigma_2, q & \stackrel{iid}{\sim} & 
        \left\{ q\frac{1}{\sqrt{2\pi}\sigma_1}
        \exp[-(y_i - \mu_1)^2/2\sigma_1^2] \right. \\
        & & ~ + ~
        \left. (1-q)\frac{1}{\sqrt{2\pi}\sigma_2}
        \exp[-(y_i - \mu_2)^2/2\sigma_2^2] \right\}
        \end{array}
        \eqno{(*)}
\eeq
We wish to estimate the parameters $\mu_1$, $\mu_2$, $\sigma_1$,
$\sigma_2$ and $q$.  Explicitly write the likelihood $g(\ul y|\mu_1,
\mu_2, \sigma_1, \sigma_2, q )$ for the observed
data $\ul y = (y_1, \ldots, y_{n})$.

\item  The likelihood is greatly simplified by using data augmentation.
  Consider the latent variable (a.k.a. missing data, or augmentation
  data)
\beq
        z_i = \left\{\begin{array}{ll}
        1 & \mbox{if $y_i$ belongs to the left hump} \\
        0 & \mbox{if $y_i$ belongs to the right hump}
        \end{array}\right. ~,~~ i=1,~\ldots,~n ~.
\eeq
% Let's say $\mu_1\leq\mu_2$ (so that $\mu_1$ is for the ``left
% hump'').  
Show that if we assume
\bea
        z_i ~|~ q & \stackrel{iid}{\sim} & Bin(1,q) \\
        y_i ~|~ \mu_1, \mu_2, \sigma_1, \sigma_2, z_i  
        & \stackrel{indep}{\sim} & 
        \left\{\begin{array}{rl}
            N(\mu_1,\sigma^2_1), & \mbox{if $z_i=1$}\\
            N(\mu_2,\sigma^2_2), & \mbox{if $z_i=0$}
            \end{array}\right.
\eea
then the distribution of the observed data $\ul y = (y_1,\ldots,
y_{n})$ is as in $(*)$.

Show, therefore, that we may take the complete-data likelihood 
for the data vectors $\ul y$ and $\ul z$ to be
\bea
        f(\ul y,\ul z|\mu_1, \mu_2, \sigma_1, \sigma_2, q)  
        & = & 
        \prod_{i=1}^n \left\{q\frac{1}{\sqrt{2\pi}\sigma_1}
        \exp[-(y_i - \mu_1)^2/2\sigma_1^2]\right\}^{z_i}  
        \\ & & ~~~ \times ~
        \left\{(1-q)\frac{1}{\sqrt{2\pi}\sigma_2}
        \exp[-(y_i - \mu_2)^2/2\sigma_2^2]\right\}^{1-z_i} 
%       & \propto &
%         \frac{q^{z_+}(1-q)^{n-z_+}}{\sigma_1^{z_+}\sigma_2^{n-z_+}}
%         \exp\{-\sum_{i=1}^{n}z_i(y_i - \mu_1)^2/2\sigma_1^2\} \\
%         & & ~ \times ~
%         \exp\{-\sum_{i=1}^{n}(1-z_i)(y_i - \mu_2)^2/2\sigma_2^2\} 
\eea
% where $n=272$ and $z_+ = \sum_{i=1}^{n} z_i$ $=$ (number of
% observations in the left hump).

\item Derive the function E-M objective function $Q(\theta,\theta^k)$
  for this problem, where $\theta=(\mu_1, \mu_2, \sigma_1, \sigma_2,
  q)$.  Are there any simplifications to the E-step?  The M-step?
  
{\bf Notes:}
\begin{itemize}
\item  $z_i ~|~ y_i, \mu_1, \mu_2, \sigma_1, \sigma_2, q$ are still
  independent Bernoulli's but now with different success
  probabilities; what are they?  
\item A modification of E-step that might be useful is to round the
  $E[z_i|y_i, \mu_1, \mu_2, \sigma_1, \sigma_2, q]$ to 0 or 1 rather
  than use the expected value itself.  How does this simplify the M-step?
\end{itemize}

\item Fit the model to the Old Faithful data, using E-M or some
  variation.  Describe your findings, focusing on the following
  questions:
  \begin{itemize}
  \item What would be a good set of starting values?
  \item How do you assess convergence?
  \item Try from several different sets of staring values.  Do you get
  different maxima?  Which maximum or maxima are ``best''?  Why?
  \end{itemize}

\item Obtain standard errors for the parameters two ways: using {\tt
  infomat.s}, and using the missing information principle.

\end{enumerate}


\end{enumerate}

\end{document}








