I'm trying to solve a homework problem and I'm not getting the answer I expect.

The problem is from Pattern Classification by Duda,Hart,Stork (problem 3.47)

Consider $\mathcal{D} = \left\{\begin{pmatrix}1\\1\end{pmatrix},\begin{pmatrix}3\\3\end{pmatrix},\begin{pmatrix}2\\*\end{pmatrix} \right\}$, sampled from a two-dimensional (separable) distribution $p(x_1,x_2)=p_{1}(x_1)p_{2}(x_2)$, with \begin{equation} p_{1}(x) \sim \begin{cases} \frac{1}{\theta_1}e^{-x/\theta_1} & \text{if}~x \geq 0\\0 & \text{otherwise}\end{cases}. \end{equation} and \begin{equation} p_{2}(x) \sim \begin{cases} \frac{1}{\theta_2} & \text{if}~0 \leq x \leq \theta_2 \\ 0 & \text{otherwise}\end{cases} \end{equation} As usual, $*$ represents a missing feature value.

(a) Start with an initial estimate of $\theta^0=\begin{pmatrix}2\\4\end{pmatrix}$ and analytically calculate $Q(\theta; \theta^0)$ - The E step of the EM algorithm. Be sure to consider the normalization of your distribution.

(b) Find the $\theta$ that maximizes your $Q(\theta; \theta^0)$ - The M step.

$Q$ is defined as

\begin{equation} Q(\theta; \theta^i) = \mathbb{E}_{\mathcal{D}_b}\left[ \ln p(\mathcal{D}_g, \mathcal{D}_b ; \theta) | \mathcal{D}_g ; \theta^i \right] \end{equation}

where $\mathcal{D}_g$ is the "good" (observed) data and $\mathcal{D}_b$ is the "bad" (missing) data. Note that $\mathcal{D}=\mathcal{D}_g \cup \mathcal{D}_b$.

I'm assuming i.i.d. sampling and we are given that the features are independent. For this reason I would expect that EM converges to the maximum likelihood estimates which is $\theta=(2,3)^T$. So I expect my solution for part (b), $\theta^1$, to be closer to $(2,3)^T$ than $\theta^0$. However this doesn't seem to be the case.


Attempt

In this problem $\mathcal{D}_g$ contains 5 observations and $\mathcal{D}_b = \{ x_{32} \}$ is the missing data. I assume the domain of $Q$ is $\theta \in \mathbb{R}^2_{++}$. Also let $\theta = (\theta_1, \theta_2)^T$ and $\theta^0 = (\theta_0^1, \theta_0^2)^T$.

\begin{align} Q(\theta; \theta^i) &= \mathbb{E}_{\mathcal{D}_b}\left[ \ln p(\mathcal{D}_g, \mathcal{D}_b ; \theta) | \mathcal{D}_g ; \theta^0 \right] \\ &= \mathbb{E}_{x_{32}}\left[ \ln p(1,1;\theta) + \ln p(3,3;\theta) + \ln p_{1}(2;\theta_1) + \ln p_{2}(x_{32};\theta_2) | \mathcal{D}_g ; \theta^0 \right] \\ &= \ln p(1,1;\theta) + \ln p(3,3;\theta) + \ln p_{1}(2;\theta_1) + \mathbb{E}_{x_{32}} \left[ \ln p_{2}(x_{32};\theta_2) | \mathcal{D}_g ; \theta^0 \right] \\ \end{align}

Consider the expression in two parts $Q(\theta; \theta^0) = f(\theta) + g(\theta; \theta^0)$

where

\begin{align} f(\theta) &= \ln p(1,1;\theta) + \ln p(3,3;\theta) + \ln p_{1}(2;\theta_1) \\ &= \begin{cases} -3 \ln \theta_1 - \frac{6}{\theta_1} - 2 \ln \theta_2 & \text{if}~\theta_2 \geq 3 \\ -\infty & \text{otherwise} \end{cases}. \end{align}

assuming $\ln(0) = -\infty$, and

\begin{align} g(\theta; \theta^0) &= \mathbb{E}_{x_{32}} \left[ \ln p_{2}(x_{32};\theta_2) | \mathcal{D}_g ; \theta^0 \right] \\ &= \int_{-\infty}^{\infty} {\left( \ln p_{2}(x_{32};\theta_2)\right) p_{2}(x_{32} | \mathcal{D}_g ; \theta^0) dx_{32}} \end{align}

Note that $x_{32}$ is independent of $\mathcal{D}_g$ because i.i.d. and $p(x_1,x_2)=p_1(x_1)p_2(x_2)$. Therefore $p_{2}(x_{32} | \mathcal{D}_g ; \theta^0) = p_{2}(x_{32} ; \theta^0_2) = 1/4$ if $x_{32} \in [0,4]$, and $0$ otherwise so then

\begin{align} g(\theta; \theta^0) &= \frac{1}{4} \int_{0}^{4} \ln p_{2}(x_{32} ; \theta_2) dx_{32} \hspace{50pt} \text{(almost surely)} \\ &= \begin{cases} -\ln \theta_2 & \text{if}~\theta_2 \geq 4 \\ -\infty & \text{otherwise} \end{cases} \end{align}

Therefore \begin{equation} Q(\theta ; \theta^0) = \begin{cases} -3\ln\theta_1 - \frac{6}{\theta_1} - 3\ln\theta_2 & \text{if}~\theta_2 \geq 4 \\ -\infty & \text{otherwise} \end{cases} \end{equation}

Clearly $\theta^1 = \text{arg}\max_{\theta} Q(\theta ; \theta^0) = (2,4)^T$.


Again this solution disagrees with my belief that $\theta^1$ should be closer to $(2,3)$ than $\theta^0$. Is this belief wrong? Have I made a mistake somewhere? I really appreciate if someone could help me understand. I think there may be a mistake in the calculation of $g(\theta ; \theta^0)$ but I can't see it.

  • Can you please explain why it has to be $\theta_2 \geq 3$ in the calculation of $f(\theta)$? – fear Jan 20 at 20:18
  • 1
    @fear Sure, it comes from the definition of $p_2$. We know that $p(3,3 ; \theta) = p_1(3; \theta_1) p_2(3; \theta_2)$. If $\theta_2 < 3$ then $p_2(3; \theta_2) = 0$, if $\theta_2 \geq 3$ then $p_2(3; \theta_2) = 1/\theta_2$. – jodag Jan 20 at 21:05
  • Thank you! One more question. Duda, in the problem solutions says that $\int_{-\infty}^{\infty} p(x_1)dx_1=1$ and thus $\int_{-\infty}^{\infty} \frac{1}{\theta_1}e^{-\frac{x_1}{\theta_1}} dx_1 =1$ and thus $\theta_1=1$. This is different from your expectation (2, 3) and also $\theta_1$ will not change values. I don't understand why $p(x_1)=\frac{1}{\theta_1}e^{-\frac{x_1}{\theta_1}}$. To my understanding it should be $p(x_1)=p_1(x_1)p_2(x_1)=\frac{1}{\theta_1}e^{-\frac{x_1}{\theta_1}} \frac{1}{\theta_2}$ and then integrate from 0 to infinity. Is there something I am missing? – fear Jan 20 at 23:12
  • @fear $\int_{-\infty}^{\infty} \frac{1}{\theta_1}e^{-\frac{x_1}{\theta_1}}dx_1 = 1$ is false. Remember that $p_1(x) = 0$ for $x < 0$ which means $1 = \int_{-\infty}^{\infty}{p_1(x_1)dx_1} = \int_{0}^{\infty}{\frac{1}{\theta_1}e^{-\frac{x_1}{\theta_1}}dx_1}$ for all $\theta_1 > 0$. – jodag Jan 21 at 0:11
  • 1
    @fear I derived a closed expression for $Q(\theta ; \theta^0)$ which is $-\infty$ for all $\theta_2 < 4$. Obviously then $\theta_2 = 3 \Rightarrow Q = -\infty$ which doesn't maximize $Q$. To maximize for $\theta_1$ just solve $\frac{\partial Q}{\partial \theta_1} = 0$ assuming that $\theta_2 \geq 4$. To maximize over $\theta_2$ note that $-3 \ln \theta_1 - \frac{6}{\theta_1} -3 \ln{\theta_2}$ is monotonically decreasing for positive $\theta_2$ so it must be $\theta_2 = 4$. – jodag Jan 21 at 0:42
up vote 3 down vote accepted

I think $g(\theta;\theta^0)$ may have a mistake. If $0 < \theta_2 < 4$, then $$g(\theta;\theta^0)=\frac{1}{4} \int_0^{\theta_2} \ln \frac{1}{\theta_2} dx_{32} = -\frac{\theta_2}{4} \ln \theta_2.$$ Thus, $Q(\theta;\theta^0) = -3 \ln \theta_1-\frac{6}{\theta_1} -(2+\frac{\theta_2}{4}) \ln \theta_2$ if $3 \le \theta_2 < 4$.

EDIT

From paper "Exercises in EM":

"E-step does not simply involve replacing missing data by their conditional expectation (although this is true for many important applications of the algorithm). Rather, the E-step takes the expected value of the complete-data log-likelihood function, conditional on the observed data. If the likelihood function takes value zero in a subset of the parameter space, then the log-likelihood function does not exist, and the EM algorithm is not applicable."

I am just replacing the missing data by the conditional expectation above and this is not correct. I think we need to be very careful with non-exponential family distributions. In the book, it mentions that "If the support does depend on $\theta$, then the monotonicity of the EM algorithm might not hold."

  • I can see that this gives me the expected solution. Can you explain why this is justified as opposed to assigning $-\infty$? – jodag Dec 11 '17 at 2:21
  • 1
    I do not see it is $-\infty$ . $g(\theta;\theta^0)=\frac{1}{4} \int_0^{4} \ln p_{2}(x_{32};\theta_2) dx_{32}$ is correct, but when $0 < \theta_2 <4$, $p_{2}(x_{32};\theta_2)$ is nonzero when $x_{32} \in (0,\theta_2)$. Thus, the integral is from 0 to $\theta_2$. – Statisfun Dec 11 '17 at 3:39
  • Maybe I'm thinking of it wrong. I'm considering $x_{32}$ as a r.v. with distribution $U(0,4)$ and $\ln p(x_{32} ; \theta_2)$ as a function of that random variable. Since $\ln p(x_{32} ; \theta_2) = -\infty$ for $x_{32} > \theta_2$ then if $\theta_2 < 4$ the function $\ln p(x_{32} ; \theta_2)$ takes on the value $-\infty$ with non-zero probability. Hense the expected value is $-\infty$ no? – jodag Dec 11 '17 at 3:43
  • Do we simply take the expected value over $\textbf{support} p(x ; \theta_2) \cap \textbf{support} p(x ; \theta_2^0)$ ? – jodag Dec 11 '17 at 3:46
  • I think you are right, sorry about my mistake. The algorithm simply got stuck at the initial values.After reading several materials, I realized that I fell into a trap. Please see my edits and thanks for asking this problem. – Statisfun Dec 11 '17 at 18:25

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.