Problem 4.16 (Maximum likelihood)

Problem 4.16 (Maximum likelihood) Assume that we are given an i.i.d. data set $(X_1,Y_1),\dots ,(X_n,Y_n)$ from the model

$$Y=f(X)+N_Y,\text{with }X\sim \mathcal{N}({\mu }_X,{\sigma }_X^2),\text{ and }{N}_Y\sim \mathcal{N}({\mu }_{N_Y},{\sigma }_{N_Y}^2)\text{ independent},$$

where the function $f$ is supposed to be known.

1. Prove that $f(x)=\mathbb{E}[Y\mid X=x]$.

2. Write $\mathbf{x}:=(x_1,\dots ,x_n)$, $\mathbf{y}:=(y_1,\dots ,y_n)$ and consider the log-likelihood function

$${\ell }_{\theta }(\mathbf{x},\mathbf{y})={\ell }_{\theta }((x_1,y_1),\dots ,(x_n,y_n))=\sum _{i=1}^n\log p_{\theta }(x_i,y_i),$$

where $p_\theta$ is the joint density over $(X,Y)$ and
$\theta := (\mu_X, \mu_{N_Y}, \sigma_X^2, \sigma_{N_Y}^2)$. Prove
that for some $c_1, c_2 \in \mathbb{R}$ with $c_2 > 0$

$$\underset{\theta }{\max }{\ell }_{\theta }(\mathbf{x},\mathbf{y})=c_2\cdot (c_1-\log \widehat{\mathrm{var}}[\mathbf{x}]-\log \widehat{\mathrm{var}}[\mathbf{y}-f(\mathbf{x})]),\text{\hspace{1em}(4.20)}$$

where $\widehat{\mathrm{var}}[\mathbf{z}]:=\frac{1}{n}{\sum }_{i=1}^n(z_i-\frac{1}{n}{\sum }_{k=1}^n{z}_k{)}^2$ estimates the variance.

Equation (4.20) motivates the comparison of expressions (4.18) and (4.19):

$$L_{X\to Y}=-\log \widehat{\mathrm{var}}[X]-\log \widehat{\mathrm{var}}[R_Y],L_{X\leftarrow Y}=-\log \widehat{\mathrm{var}}[R_X]-\log \widehat{\mathrm{var}}[Y],$$

with $R_Y:=Y-{\hat{f}}_Y(X)$ and $R_X:=X-{\hat{f}}_X(Y)$ the regression residuals. The main difference is that in this exercise, we have used the conditional mean and not the outcome of the regression method.

Solution a) As stated this needs ${\mu }_{N_Y}=0$: by linearity of conditional expectation,

$$\mathbb{E}[Y\mid X=x]=\mathbb{E}[f(X)\mid X=x]+\mathbb{E}[N_Y\mid X=x]=f(x)+{\mu }_{N_Y}.$$

Why? Let us use Part A of cond_exp.pdf. Set ${\kappa }_x:={\delta }_{f(x)}$. This is a regular conditional distribution of $f(X)$ given $X$: (K1), (K2) are immediate, and (K3) holds since for Borel $A,B$

$${\int }_A{\delta }_{f(x)}(B)P_X(dx)={\int }_A{1}_{f^{-1}(B)}(x)P_X(dx)=P_X(A\cap f^{-1}(B))=P(X\in A,f(X)\in B).$$

Hence we can write

$$E\left[f(X)\mid X=x\right]=\int z{\delta }_{f(x)}(dz)=f(x).$$

And the second equals the unconditional mean $\mathbb{E}[N_Y]$ because $N_Y\perp \perp X$, because of lemma 2 (note that the independence is used in second to last equation in the equation chain). So $f(x)=\mathbb{E}[Y\mid X=x]$ iff ${\mu }_{N_Y}=0$. This is a harmless normalization: replacing $(f,{\mu }_{N_Y})$ by $(\tilde{f},0)$ with $\tilde{f}:=f+{\mu }_{N_Y}$ gives the identical joint distribution, so only the sum $f+{\mu }_{N_Y}$ — that is, the regression function $x\mapsto \mathbb{E}[Y\mid X=x]$ — is identifiable.

Solution b) The joint density factorizes as $p_{\theta }(x,y)=p(x)p(y\mid x)$ with $X\sim \mathcal{N}({\mu }_X,{\sigma }_X^2)$ and $Y\mid X=x\sim \mathcal{N}(f(x)+{\mu }_{N_Y},{\sigma }_{N_Y}^2)$.

Why does the conditional distribution look like that? Left as a gap.

Writing $r_k:=y_k-f(x_k)$,

$${\ell }_{\theta }(\mathbf{x},\mathbf{y})=\sum _{k=1}^n\left[-\frac{1}{2}\log (2\pi {\sigma }_X^2)-\frac{(x_k-{\mu }_X{)}^2}{2{\sigma }_X^2}\right]+\sum _{k=1}^n\left[-\frac{1}{2}\log (2\pi {\sigma }_{N_Y}^2)-\frac{(r_k-{\mu }_{N_Y}{)}^2}{2{\sigma }_{N_Y}^2}\right]=:{\ell }_{\theta }(\mathbf{x})+{\ell }_{\theta }(\mathbf{r}).$$

Note the $-\frac{n}{2}\log (2\pi {\sigma }^2)$ terms are not constant in $\theta$; without them the ${\sigma }^2$-derivatives would be strictly positive and have no critical point. Note also that ${\mu }_{N_Y}$ must stay free: $\theta$ has four components, and maximizing over the slice ${\mu }_{N_Y}=0$ is a different (smaller) number in general.

The two blocks have the same form, so it suffices to maximize one, say ${\ell }_{\theta }(\mathbf{x})$ over $({\mu }_X,{\sigma }_X^2)$. First-order conditions:

$$\frac{\partial \ell }{\partial {\mu }_X}=\sum _k\frac{x_k-{\mu }_X}{{\sigma }_X^2}=0\Rightarrow {\hat{\mu }}_X=\bar{x}:=\frac{1}{n}{\sum }_k{x}_k,$$ $$\frac{\partial \ell }{\partial {\sigma }_X^2}=-\frac{n}{2{\sigma }_X^2}+\sum _k\frac{(x_k-{\mu }_X{)}^2}{2{\sigma }_X^4}=0\Rightarrow {\hat{\sigma }}_X^2=\frac{1}{n}{\sum }_k(x_k-\bar{x}{)}^2=\widehat{\mathrm{var}}[\mathbf{x}].$$

(This critical point is the maximum, because the density funciton is smooth and can be checked by having the Hessian minors’ determinants have signs -1, 1. Then we can check that it is a global min) Substituting back, the quadratic term collapses to a constant:

$$\sum _k\frac{(x_k-\bar{x}{)}^2}{2\widehat{\mathrm{var}}[\mathbf{x}]}=\frac{n\widehat{\mathrm{var}}[\mathbf{x}]}{2\widehat{\mathrm{var}}[\mathbf{x}]}=\frac{n}{2},$$

so

$$\underset{{\mu }_X,{\sigma }_X^2}{\max }{\ell }_{\theta }(\mathbf{x})=-\frac{n}{2}\log \widehat{\mathrm{var}}[\mathbf{x}]-\frac{n}{2}\log (2\pi )-\frac{n}{2}.$$

The residual block is identical with $r_k$ in place of $x_k$: ${\hat{\mu }}_{N_Y}=\bar{r}$, ${\hat{\sigma }}_{N_Y}^2=\widehat{\mathrm{var}}[\mathbf{y}-f(\mathbf{x})]$ — the free ${\mu }_{N_Y}$ is exactly what centers the residuals, producing the centered $\widehat{\mathrm{var}}$ rather than $\frac{1}{n}{\sum }_k{r}_k^2$: the ${\sigma }^2$-condition is ${\hat{\sigma }}_{N_Y}^2=\frac{1}{n}{\sum }_k(r_k-{\mu }_{N_Y}{)}^2$; with ${\mu }_{N_Y}$ free, ${\hat{\mu }}_{N_Y}=\bar{r}$ and this is $\widehat{\mathrm{var}}[\mathbf{r}]$; with ${\mu }_{N_Y}$ fixed at ${\mu }_0$ it would be $\widehat{\mathrm{var}}[\mathbf{r}]+(\bar{r}-{\mu }_0{)}^2$ and (4.20) would fail as stated. Adding and factoring out $n/2$:

$$\underset{\theta }{\max }{\ell }_{\theta }(\mathbf{x},\mathbf{y})=\frac{n}{2}(\underset{c_1}{\underbrace{-2\log (2\pi )-2}}-\log \widehat{\mathrm{var}}[\mathbf{x}]-\log \widehat{\mathrm{var}}[\mathbf{y}-f(\mathbf{x})]),$$

which is (4.20) with $c_2=n/2>0$ and $c_1=-2(\log 2\pi +1)$. The coefficient 1 on each log term pins $c_2=n/2$ uniquely — $(c_1,c_2)$ cannot be rescaled to, say, $(n,-\log 2\pi -1)$ (why?). Since $c_1,c_2$ do not depend on the direction of the model, comparing causal directions by maximum likelihood reduces to comparing $\log \widehat{\mathrm{var}}[\mathbf{x}]+\log \widehat{\mathrm{var}}[\mathbf{y}-f(\mathbf{x})]$ with its flipped counterpart, i.e. the scores (4.18)/(4.19).