Math Problems Vault

Problem 6.4.1

Let

u(t)= \begin{cases} c \int _{0}^{t} e^{-1 / x(1-x)} d x & \text{ for } 0<t<1 \\ 0 & \text{ for } t \leq 0 \\ 1 & \text{ for } t \geq 1\end{cases}

Show that for a suitable $c$ , the function $u$ is continuous and infinitely differentiable for $-\infty <t<\infty$ .

?

Problem 6.4.2

Show that the density (4.1) with $\Omega =(0, \infty )$ satisfies all conditions of Theorem 3.10 with the exception of (d) of Theorem 2.6.

?

Problem 6.4.3

Show that the density (4.4) with $\Omega =(0, \infty )$ satisfies all conditions of Theorem 3.10.

?

Problem 6.4.4

In Example 4.5, evaluate the estimators (4.8) and (4.14) for the Cauchy case, using for $\tilde{\theta }_{n}$ the sample median.

?

Problem 6.4.5

In Example 4.7, show that $l(\theta )$ is concave.

?

Problem 6.4.6

In Example 4.7, if $\eta =\xi$ , show how to obtain a $\sqrt{n}$ -consistent estimator by equating sample and population second moments.

?

Problem 6.4.7

In Theorem 4.8, show that $\sigma_{11}=\sigma_{12}$ .

?

Problem 6.4.8

Without using Theorem 4.8, in Example 4.13 show that the EM sequence converges to the MLE.

?

Problem 6.4.9

(a)(b)(c)

Consider the following 12 observations from a bivariate normal distribution with parameters $\mu_{1}=\mu_{2}=0, \sigma_{1}^{2}, \sigma_{2}^{2}, \rho$ :

$x_{1}$	1	1	-1	-1	2	2	-2	-2	*	*	*	*
$x_{2}$	1	-1	1	-1	*	*	*	*	2	2	-2	-2

where " $*$ " represents a missing value.

?

(a)

Show that the likelihood function has global maxima at $\rho = \pm 1 / 2, \sigma_{1}^{2}=\sigma_{2}^{2}=8 / 3$ , and a saddlepoint at $\rho =0, \sigma_{1}^{2}=\sigma_{2}^{2}=5 / 2$ .

(b)

Show that if an EM sequence starts with $\rho =0$ , then it remains with $\rho =0$ for all subsequent iterations.

(c)

Show that if an EM sequence starts with $\rho$ bounded away from zero, it will converge to a maximum.

Remark.

?

This problem is due to Murray (1977), and is discussed by Wu (1983).

Problem 6.4.10

Show that if the EM complete-data density $f(\mathbf{y}, \mathbf{z} \mid \theta )$ of (4.21) is in a curved exponential family, then the hypotheses of Theorem 4.12 are satisfied.

?

Problem 6.4.11

(a)(b)

In the EM algorithm, calculation of the E-step, the expectation calculation, can be complicated. In such cases, it may be possible to replace the E-step by a Monte Carlo evaluation, creating the MCEM algorithm (Wei and Tanner 1990). Consider the following MCEM evaluation of $Q\left(\theta \mid \hat{\theta }_{(j)}, \mathbf{y}\right)$ : Given $\hat{\theta }_{(j)}^{(k)}$

Generate $Z_{1}, \cdots , Z_{k}$ , iid, from $k\left(\mathbf{z} \mid \hat{\theta }_{(j)}^{(k)}, \mathbf{y}\right)$ ,
Let $\hat{Q}\left(\theta \mid \hat{\theta }_{(j)}^{(k)}, \mathbf{y}\right)=\frac{1}{k} \sum_{i=1}^{k} \log L(\theta \mid \mathbf{y}, \mathbf{z})$

and then calculate $\hat{\theta }_{(j+1)}^{(k)}$ as the value that maximizes $\hat{Q}\left(\theta \mid \hat{\theta }_{(j)}^{(k)}, \mathbf{y}\right)$ .

?

(a)

Show that $\hat{Q}\left(\theta \mid \hat{\theta }_{(j)}^{(k)}, \mathbf{y}\right) \rightarrow \hat{Q}\left(\theta \mid \hat{\theta }_{(j)}, \mathbf{y}\right)$ as $k \rightarrow \infty$ .

(b)

What conditions will ensure that $L\left(\hat{\theta }_{(j+1)}^{(k)} \mid \mathbf{y}\right) \geq L\left(\hat{\theta }_{(j)}^{(k)} \mid \mathbf{y}\right)$ for sufficiently large $k$ ? Are the hypotheses of Theorem 4.12 sufficient?

Problem 6.4.12

(a)(b)(c)

For the mixture distribution of Example 4.7, that is,

X_{i} \sim \theta g(x)+(1-\theta ) h(x), \quad i=1, \ldots , n, \text{ independent }

where $g(\cdot )$ and $h(\cdot )$ are known, an EM algorithm can be used to find the ML estimator of $\theta$ . Let $Z_{1}, \cdots , Z_{n}$ , where $Z_{i}$ indicates from which distribution $X_{i}$ has been drawn, so

\begin{align} & X_{i} \mid Z_{i}=1 \sim g(x) \\ & X_{i} \mid Z_{i}=0 \sim h(x). \end{align}

?

(a)

Show that the complete-data likelihood can be written

L(\theta \mid \mathbf{x}, \mathbf{z})=\prod _{i=1}^{n}\left[z_{i} g\left(x_{i}\right)+\left(1-z_{i}\right) h\left(x_{i}\right)\right] \theta ^{z_{i}}(1-\theta )^{1-z_{i}}.

(b)

Show that $\mathbb {E}\left[Z_{i} \mid \theta , x_{i}\right]=\theta g\left(x_{i}\right) /\left[\theta g\left(x_{i}\right)+(1-\theta ) h\left(x_{i}\right)\right]$ and hence that the EM sequence is given by

\hat{\theta }_{(j+1)}=\frac{1}{n} \sum _{i=1}^{n} \frac{\hat{\theta }_{(j)} g\left(x_{i}\right)}{\hat{\theta }_{(j)} g\left(x_{i}\right)+\left(1-\hat{\theta }_{(j)}\right) h\left(x_{i}\right)}.

(c)

Show that $\hat{\theta }_{(j)} \rightarrow \hat{\theta }$ , the ML estimator of $\theta$ .

Problem 6.4.13

For the situation of Example 4.10:

?

(a)

Show that the M-step of the EM algorithm is given by

\hat{\mu }=\left(\sum _{i=1}^{4} \sum _{j=1}^{n_{i}} y_{i j}+z_{1}+z_{2}\right) / 12,

\begin{align} \hat{\alpha }_{i} & =\left(\sum _{j=1}^{2} y_{i j}+z_{i}\right) / 3-\hat{\mu }, \quad i=1,3 \\ & =\left(\sum _{j=1}^{3} y_{i j}\right) / 3-\hat{\mu }, \quad i=2,4. \end{align}

(b)

Show that the E-step of the EM algorithm is given by

z_{i}=\mathbb {E}\left[Y_{i 3} \mid \mu =\hat{\mu }, \alpha _{i}=\hat{\alpha }_{i}\right]=\hat{\mu }+\hat{\alpha }_{i} \quad i=1,3.

(c)

Under the restriction $\sum_{i} \alpha_{i}=0$ , show that the EM sequence converges to $\hat{\alpha }_{i}= \bar{y}_{i.}-\hat{\mu }$ , where $\hat{\mu }=\sum_{i} \bar{y}_{i.} / 4$ .

(d)

Under the restriction $\sum_{i} n_{i} \alpha_{i}=0$ , show that the EM sequence converges to $\hat{\alpha }_{i}= \bar{y}_{i.}-\hat{\mu }$ , where $\hat{\mu }=\sum_{i j} y_{i j} / 10$ .

(e)

For a general one-way layout with $a$ treatments and $n_{i j}$ observations per treatment, show how to use the EM algorithm to augment the data so that each treatment has $n$ observation. Write down the EM sequence, and show what it converges to under the restrictions of parts (c) and (d).

[The restrictions of parts (c) and (d) were encountered in Example 3.4.9, where they led, respectively, to an unweighted means analysis and a weighted means analysis.]

Problem 6.4.14

(a)(b)(c)

In the two-way layout (see Example 3.4.11), the EM algorithm can be very helpful in computing ML estimators in the unbalanced case. Suppose that we observe

Y_{i j k}: N\left(\xi _{i j}, \sigma ^{2}\right), \quad i=1, \ldots , I, \quad j=1, \ldots , J, \quad k=1, \ldots , n_{i j},

where $\xi_{i j}=\mu +\alpha_{i}+\beta_{j}+\gamma_{i j}$ . The data will be augmented so that the complete data have $n$ observations per cell.

?

(a)

Show how to compute both the E-step and the M-step of the EM algorithm.

(b)

Under the restriction $\sum_{i} \alpha_{i}=\sum_{j} \beta_{j}=\sum_{i} \gamma_{i j}=\sum_{j} \gamma_{i j}=0$ , show that the EM sequence converges to the ML estimators corresponding to an unweighted means analysis.

(c)

Under the restriction $\sum_{i} n_{i}. \alpha_{i}=\sum_{j} n_{\cdot j} \beta_{j}=\sum_{i} n_{i \cdot } \gamma_{i j}=\sum_{j} \cdot j \gamma_{i j}=0$ , show that the EM sequence converges to the ML estimators corresponding to a weighted means analysis.

Problem 6.4.15

For the one-way layout with random effects (Example 3.5.1), the EM algorithm is useful for computing ML estimates. (In fact, it is very useful in many mixed models; see Searle et al. 1992, Chapter 8.) Suppose we have the model

X_{i j}=\mu +A_{i}+U_{i j} \quad \left(j=1, \ldots , n_{i}, i=1, \ldots , s\right)

where $A_{i}$ and $U_{i j}$ are independent normal random variables with mean zero and known variance. To compute the ML estimates of $\mu , \sigma_{U}^{2}$ , and $\sigma_{U}^{2}$ it is typical to employ an EM algorithm using the unobservable $A_{i}$ 's as the augmented data. Write out both the E-step and the M-step, and show that the EM sequence converges to the ML estimators.

?

Problem 6.4.16

Maximum likelihood estimation in the probit model of Section 3.6 can be implemented using the EM algorithm. We observe independent Bernoulli variables $X_{1}, \ldots , X_{n}$ , which depend on unobservable variables $Z_{i}$ distributed independently as $N\left(\zeta_{i}, \sigma^{2}\right)$ , where

X_{i}= \begin{cases} 0 & \text{ if } Z_{i} \leq u \\ 1 & \text{ if } Z_{i}>u.\end{cases}

Assuming that $u$ is known, we are interested in obtaining ML estimates of $\zeta$ and $\sigma^{2}$ .

?

(a)

Show that the likelihood function is $p^{\Sigma x_{i}}(1-p)^{n-\Sigma x_{i}}$ , where

p=\mathbb {P}\left(Z_{i}>u\right)=\Phi \left(\frac{\zeta -u}{\sigma }\right).

(b)

If we consider $Z_{1}, \ldots , Z_{n}$ to be the complete data, the complete-data likelihood is

\prod _{i=1}^{n} \frac{1}{\sqrt{2 \pi } \sigma } e^{-\frac{1}{2 \sigma ^{2}}\left(z_{i}-\zeta \right)^{2}}

and the expected complete-data log likelihood is

-\frac{n}{2} \log \left(2 \pi \sigma ^{2}\right)-\frac{1}{2 \sigma ^{2}} \sum _{i=1}^{n}\left[\mathbb {E}\left[Z_{i}^{2} \mid X_{i}\right]-2 \zeta \mathbb {E}\left[Z_{i} \mid X_{i}\right]+\zeta ^{2}\right].

(c)

Show that the EM sequence is given by

\begin{align} & \hat{\zeta }_{(j+1)}=\frac{1}{n} \sum _{i=1}^{n} t_{i}\left(\hat{\zeta }_{(j)}, \hat{\sigma }_{(j)}^{2}\right) \\ & \hat{\sigma }_{(j+1)}^{2}=\frac{1}{n}\left[\sum _{i=1}^{n} v_{i}\left(\hat{\zeta }_{(j)}, \hat{\sigma }_{(j)}^{2}\right)-\frac{1}{n}\left(\sum _{i=1}^{n} t_{i}\left(\hat{\zeta }_{(j)}, \hat{\sigma }_{(j)}^{2}\right)\right)^{2}\right] \end{align}

where

t_{i}\left(\zeta , \sigma ^{2}\right)=\mathbb {E}\left[Z_{i} \mid X_{i}, \zeta , \sigma ^{2}\right] \quad \text{ and } \quad v_{i}\left(\zeta , \sigma ^{2}\right)=\mathbb {E}\left[Z_{i}^{2} \mid X_{i}, \zeta , \sigma ^{2}\right].

(d)

Show that

\begin{align} & \mathbb {E}\left[Z_{i} \mid X_{i}, \zeta , \sigma ^{2}\right]=\zeta +\sigma H_{i}\left(\frac{u-\zeta }{\sigma }\right) \\ & \mathbb {E}\left[Z_{i}^{2} \mid X_{i}, \zeta , \sigma ^{2}\right]=\zeta ^{2}+\sigma ^{2}+\sigma (u+\zeta ) H_{i}\left(\frac{u-\zeta }{\sigma }\right) \end{align}

where

H_{i}(t)= \begin{cases} \frac{\varphi (t)}{1-\Phi (t)} & \text{ if } X_{i}=1 \\ -\frac{\varphi (t)}{\Phi (t)} & \text{ if } X_{i}=0\end{cases}

(e)

Show that $\hat{\zeta }_{(j)} \rightarrow \hat{\zeta }$ and $\hat{\sigma }_{(j)}^{2} \rightarrow \hat{\sigma }^{2}$ , the ML estimates of $\zeta$ and $\sigma^{2}$ .

Problem 6.4.17

Verify (4.30).

?

Problem 6.4.18

(a)(b)(c)

The EM algorithm can also be implemented in a Bayesian hierarchical model to find a posterior mode. Recall the model (4.5.5.1),

\begin{align} & X \mid \theta \sim f(x \mid \theta ), \\ & \Theta \mid \lambda \sim \pi (\theta \mid \lambda ), \\ & \Lambda \sim \gamma (\lambda ), \end{align}

where interest would be in estimating quantities from $\pi (\theta \mid x)$ . Since

\pi (\theta \mid x)=\int \pi (\theta , \lambda \mid x) d \lambda ,

where $\pi (\theta , \lambda \mid x)=\pi (\theta \mid \lambda , x) \pi (\lambda \mid x)$ , the EM algorithm is a candidate method for finding the mode of $\pi (\theta \mid x)$ , where $\lambda$ would be used as the augmented data.

?

(a)

Define $k(\lambda \mid \theta , x)=\pi (\theta , \lambda \mid x) / \pi (\theta \mid x)$ , and show that

\log \pi (\theta \mid x)=\int \log \pi (\theta , \lambda \mid x) k\left(\lambda \mid \theta ^{*}, x\right) d \lambda -\int \log k(\lambda \mid \theta , x) k\left(\lambda \mid \theta ^{*}, x\right) d \lambda .

(b)

If the sequence $\left\{ \hat{\theta }_{(j)}\right\}$ satisfies

\max _{\theta } \int \log \pi (\theta , \lambda \mid x) k\left(\lambda \mid \theta _{(j)}, x\right) d \lambda =\int \log \pi \left(\theta _{(j+1)}, \lambda \mid x\right) k\left(\lambda \mid \theta _{(j)}, x\right) d \lambda ,

show that $\log \pi \left(\theta_{(j+1)} \mid x\right) \geq \log \pi \left(\theta_{(j)} \mid x\right)$ . Under what conditions will the sequence $\left\{ \hat{\theta }_{(j)}\right\}$ converge to the mode of $\pi (\theta \mid x)$ ?

(c)

For the hierarchy

\begin{align} & X \mid \theta \sim N(\theta ), 1), \\ & \Theta \mid \lambda \sim N(\lambda , 1)), \\ & \Lambda \sim \operatorname {Uniform}(-\infty , \infty ), \end{align}

show how to use the EM algorithm to calculate the posterior mode of $\pi (\theta \mid x)$ .

Problem 6.4.19

(a)(b)

There is a connection between the EM algorithm and Gibbs sampling, in that both have their basis in Markov chain theory. One way of seeing this is to show that the incomplete-data likelihood is a solution to the integral equation of successive substitution sampling (see Problems 4.5.9-4.5.11), and that Gibbs sampling can then be used to calculate the likelihood function. If $L(\theta \mid \mathbf{y})$ is the incomplete-data likelihood and $L(\theta \mid \mathbf{y}, \mathbf{z})$ is the complete-data likelihood, define

\begin{align} L^{*}(\theta \mid \mathbf{y}) & =\frac{L(\theta \mid \mathbf{y})}{\int L(\theta \mid \mathbf{y}) d \theta } \\ L^{*}(\theta \mid \mathbf{y}, \mathbf{z}) & =\frac{L(\theta \mid \mathbf{y}, \mathbf{z})}{\int L(\theta \mid \mathbf{y}, \mathbf{z}) d \theta } \end{align}

?

(a)

Show that $L^{*}(\theta \mid \mathbf{y})$ is the solution to

L^{*}(\theta \mid \mathbf{y})=\int \left[\int L^{*}(\theta \mid \mathbf{y}, \mathbf{z}) k\left(\mathbf{z} \mid \theta ^{\prime }, \mathbf{y}\right) d \mathbf{z}\right] L^{*}\left(\theta ^{\prime } \mid \mathbf{y}\right) d \theta ^{\prime }

where, as usual, $k(\mathbf{z} \mid \theta , \mathbf{y})=L(\theta \mid \mathbf{y}, \mathbf{z}) / L(\theta \mid \mathbf{y})$ .

(b)

Show how the sequence $\theta_{(j)}$ from the Gibbs iteration,

\begin{align} & \theta _{(j)} \sim L^{*}\left(\theta \mid \mathbf{y}, \mathbf{z}_{(j-1)}\right) \\ & \mathbf{z}_{(j)} \sim k\left(\mathbf{z} \mid \theta _{(j)}, \mathbf{y}\right) \end{align}

will converge to a random variable with density $L^{*}(\theta \mid \mathbf{y})$ as $j \rightarrow \infty$ . How can this be used to compute the likelihood function $L(\theta \mid \mathbf{y})$ ?

[Using the functions $L(\theta \mid \mathbf{y}, \mathbf{z})$ and $k(\mathbf{z} \mid \theta , \mathbf{y})$ , the EM algorithm will get us the ML estimator from $L(\theta \mid \mathbf{y})$ , whereas the Gibbs sampler will get us the entire function. This likelihood implementation of the Gibbs sampler was used by Casella and Berger (1994) and is also described by Smith and Roberts (1993). A version of the EM algorithm, where the Markov chain connection is quite apparent, was given by Baum and Petrie (1966) and Baum et al. (1970).]