\documentclass[12pt]{article}
\input{defs}

\usepackage{amsmath,color,graphicx}

\setlength{\parindent}{0pt}
\setlength{\parskip}{4pt}

\addtolength{\evensidemargin}{-0.5in}
\addtolength{\oddsidemargin}{-0.5in}
\addtolength{\textwidth}{1in}
\addtolength{\topmargin}{-0.5in}
\addtolength{\textheight}{1in}

\newcommand{\mybox}{\fbox{\rule[-4pt]{0pt}{16pt}\hspace{24pt}}}

\renewcommand{\today}{Thursday, September 19}
\title{\Large\bf 6.867 Machine Learning}
\author{\bf Problem set 1 --- Solutions}

\begin{document}
\maketitle

\input{../guide.tex}

Please address questions and comments to \texttt{6.867-staff@ai.mit.edu}.

\section*{Problem 1: regression}

{\bf Reference:} Lecture two; Chapter five.


Here we will be using a regression method to predict 
mileage per gallon (MPG) of gas for different vehicles. You'll find the data in the file
``mpg.dat''. Information about the data, including the column
interpretation can be found in the file ``mpg.names''. These
files, like many other data files in the course, are taken from the
UCI Machine Learning Repository
\texttt{http://www.ics.uci.edu/$\sim$mlearn/MLSummary.html}.


\begin{enumerate}

\item[] First, we will use a linear regression model to predict
the MPG values, using squared-error as the criterion to minimize.  In
other words $y=f(\mathbf{x};\hat{\mathbf{w}})=\hat{w}_0 +
\sum_{i=1}^{7} \hat{w}_i x_{i}$, where $\hat{\mathbf{w}} = \arg
\min_{\textbf{w}} \sum_{t=1}^{n} (\mathbf{y_t} -
f(\mathbf{x_t};\mathbf{w}))^2$; here $y_t$ is the MPG of a car $t$,
$x_t$ is the 7-component input vector representing the car $t$, and
$n$ is the number of training examples.

To allow experiments with linear regression, the following Matlab
functions have been provided.
\begin{verbatim}
w = train(trainX,trainY)
\end{verbatim} 

given the set of training
vectors ${\bf x}_t,\:t=1,\ldots,n$ (rows of \texttt{trainX}) and
corresponding output values $y_t,\:t=1,\ldots,n$ 
(elements of \texttt{trainY}),  returns in \texttt{w} 
the ``optimal'' weight vector $\hat{\mathbf{w}}$. 

\begin{verbatim}
Y = predict(X,w)
\end{verbatim}

given a set of input vectors (rows of \texttt{X}) and a weight vector
$\hat{\mathbf{w}}$, returns in elements of \texttt{Y}
output values $\hat{y}$ predicted for \texttt{X} using $\hat{\mathbf{w}}$.

\begin{verbatim}
[trainE, testE] = trainNtest(trainX,trainY,testX,testY)
\end{verbatim}

given the input and output values in a training set and a test set,
trains and evaluates linear regression model, and returns in \texttt{trainE} and \texttt{testE} the training and test
errors, respectively.

\begin{verbatim}
[trainE,testE] = testLinear(x,y,numtrain)
\end{verbatim}

takes as input the complete data set (both input vectors and output values), and the desired training set size, and
returns the mean training and test errors.

To test our linear regression model, we will partition the available
data into disjoint training and test sets. This will be done, for a
training set size $N$, by using the first $N$ lines in the data file
for training, and the
remaining lines for testing.

The easiest way to read in the data in Matlab, and 
\begin{verbatim}
data = load('mpg.dat');
x = data(:,1:7);
y = data(:,8);
\end{verbatim}

\begin{figure}
\centering{
\includegraphics[height=.4\textheight]{plot1}
\caption{Training / test error of linear regression as a function of
training set size}\label{fig:plot1}
}
\end{figure}


The training and test errors of linear regression are plotted in
Figure \ref{fig:plot1} as a function of an increasing size of
the training set. Make sure you get the same values.  

\problem{8} How would you expect the graphs on Figure \ref{fig:plot1}
to behave with an unlimited supply of examples? In this case we
estimate the regression model with only a small training set as before
but the size of the test set would be infinite. In your answer refer
to the values of the training and test errors as well as to
their behavior with respect to each other.

\begin{answer}
Note that both the training error and the test error can be viewed as
estimates of the true generalization error. Both are consistent
estimates; as the size of the data increases, their value approaches
generalization error. In this case, the test error will approach that
limit. Its value is likely to be somewhat larger than the training
error, since the training error on data set of fixed (small) size is
``optimistic'' - the model parameters are optimized to minimize the
training error. If the gap between the test error and the training
error is large we say that the learning algorithm is overfitting.
\end{answer}

\problem{5} In all but very unlikely scenarios the value of the
training error on a training set of size 5 in MPG should not depend on
the actual data. Explain why, and write down that value. Do you get
this answer by running \texttt{testLinear} on the training set of size
5? If no, explain why.

\begin{answer}
The linear regression model for MPG has 8 parameters
$w_0,\ldots,w_7$. Fitting the model to 5 data points can be expressed
as the system of 5 linear equations and 8 unknowns: 
\[\left\{\begin{array}{l}
y1=w_0+w_1x_{1,1}+\ldots+w_7x_{1,7}\\
\vdots\\
y5=w_0+w_1x_{5,1}+\ldots+w_7x_{5,7}
\end{array}\right.
\]
where $x_{i,j}$ denotes the $j$-th component of the $i$-th input
vector. This system will usually have an exact solution - in fact,
if the points are in general position (not collinear), there will be an
infinite number of solutions. The training error corresponding to such
a solution is zero.

The ``unlikely scenario'' to which we refer in the question is the
scenario in which the MPG measurements taken for two (or more) cars
with identical parameters produce different values. In that case, the
system above would not be consistent, and parameters of the linear
regression model would result in a non-zero training error. As we get
more training points the system of equations would typically be 
``inconsistent'' (more equations than unknowns). 

In practice, finding the parameters involves a number of calculation
steps, such as finding a matrix pseudo-inverse, that inevitably
introduce numerical error (due to finite precision). Therefore, our
Matlab program reports a very small, but non-zero training error; we
got an error of the order $1o^{-28}$.
\end{answer}

\item[] We will now move on to polynomial regression. We will predict
the MPG values using a function of the form:
$$f(\mathbf{x};\mathbf{w}) = w_0 + \sum_{i=1}^{7} \sum_{d=1}^m
w_{i,d} x_i^d,$$
where again, the weights $\mathbf{w}$ are chosen so as to minimize
the mean squared error on the training set. Think about why we also
include all lower order polynomial terms up to the highest order
rather than just the highest ones [do not turn in an answer]. 

Note that we only use features which are powers of a single input
feature. We do so mostly in order to simplify the problem. In most
cases, it is more beneficial to use features which are products of
different input features, and perhaps also their powers. Think of why
such features are usually more powerful [do not turn in an
answer].

You can use the provided Matlab function
\begin{verbatim}
xx = degexpand(x,m)
\end{verbatim}
that for given set of input vectors \texttt{x} and a degree \texttt{m}
returns the set of
expanded vectors (containing all the powers of the original inputs).

You will have to slightly modify \texttt{testLinear} in order to 
experiment with polynomial regression. You do not have to turn in the
modified code.

{\bf NOTE:} When the degree is high, some of the features will have
extremely high values, while others will have very low values. This
causes severe numeric precision problems with matrix inversion, and
yields wrong answers. To overcome this problem, one must 
appropriately scale each feature $x_i^d$ included in the regression
model, to bring all features to roughly the same magnitude. Our
version
of \texttt{degexpand}, for instance, normalizes each feature by the maximum absolute value of that feature.

Of course, this scaling of features can be safely done only if
it does not change the regression predictions.  That is, given
training feature vectors and output values
$\{\mathbf{x}_t,y_t\}_{t=1,\ldots,n}$ and a test input vector
$\mathbf{x}_{\text{test}}$, and scaling factors $\{ \alpha_i
\}_{i=1,\ldots,7}$, we would like to make sure that the prediction of
the test output value would be the same if we trained a linear
regression on $\{\mathbf{\tilde{x}}_t,y_t\}_{t=1,\ldots,n}$, where
$\tilde{x}_{t,i} = \alpha_i x_{t,i}$, and predicted on
$\mathbf{\tilde{x}}_{\text{test}}$, where $\tilde{x}_{\text{test},i} =
\alpha_i x_{\text{test},i}$. In fact, it is enough to prove this for
the linear model (maximum degree one). The result extends to scaling
each ``extended'' feature $x_i^d$ (which is what we actually do),
since this is just linear regression using these ``extended''
features.  Now, let's assume $\hat{w}_i$ is the ``optimal'' weight
corresponding to the $i^{th}$ feature when we use the original
non-scaled inputs. Then, $\tilde{w}_i = \hat{w}_i/\alpha_i$ gives the
``optimal'' weights for the scaled features since $\sum_{i=1}^7
\tilde{w}_i \tilde{x}_i = \sum_{i=1}^7 \hat{w}_i x_i$. In other words,
we make the same predictions in both cases.  More generally, since
$\mathbf{\tilde{x}}$ is an (invertible) linear function of
$\mathbf{x}$, the set $A$ of linear function of $\mathbf{x}$ is
exactly equal to the set $\tilde{A}$ of linear functions of
$\mathbf{\tilde{x}}$. Therefore, the optimal predictor for the
original problem, $\argmin_{w\in A}(\text{training error})$ is exactly
the same as the optimal predictor $\argmin_{w\in
\tilde{A}}(\text{training error})$ for the scaled problem.

\begin{advice}
In order to understand the expansion that \texttt{degexpand}, apply it
on a small data set, for example, on
\begin{verbatim}
x=[1,   2,   3;
   2,   0.5, 1.5;
   0,   1,   2];
\end{verbatim}
and make sure you understand the results of \texttt{degexpnad(x)}.
\end{advice}


\problem{7} For a training set size of 300, turn in the mean squared
training and test errors for maximal degrees of zero through eight.

``Sanity check'': for maximal degree two, you should get a training error of
$6.8164$ and a test error of $9.4656$).

\begin{answer}

The training and test errors we got are listed in the table below:
\centering
\begin{small}
\begin{tabular}{l|c|c}
\hline
Degree & Training error & Test error \\
\hline
0 & 63.887201 & 50.599153 \\
1 & 10.222106 & 13.702433 \\
2 & 6.816417 & 9.465568 \\
3 & 6.383646 & \textbf{7.042212} \\
4 & 6.200593 & 11.294545 \\
5 & 6.055586 & 105.603196 \\
6 & 5.973245 & 13714127.434883 \\
7 & 5.598525 & 418342920955.738220 \\
8 & 5.464701 & 254846135855.677948 \\
\hline
\end{tabular}
\end{small}
\end{answer}


\problem{15} Explain the qualitative behavior of the test error as a function
of the polynomial degree. Which degree seems to be the best choice?

\begin{answer}
The higher the maximal degree in the polynomial regression model, the
better it can fit complex patterns in the data. The model classes
indexed by the maximal degree are also nested in the sense that the
models with higher maximal degree include those with lower maximal
degree as special cases (we can simply set the higher order
coefficients to zero). This means that by increasing the maximal
degree, we never loose any ability to fit the data, only increase it.
This accounts for the monotonically decreasing value of the training
error as we increase the degree. Initially, this also leads to an
improved generalization performance (as measured by the test error) -
the model fits better the true underlying distribution of the
problem. However, the more complex model also has more parameters, and
consequently requires more data points to adequately estimate those
parameters. Without providing additional data, we are bound to start
overfitting the complex model to the data - that is, instead of
fitting the true behavior, the model fits ``artifacts'' that are
present in the specific data set but are not typical for the data in
general. As a result, the generalization error (and its estimate, the
test error) start deteriorating.

The best choice of a degree here is to pick the degree with minimal
test error, which is  3. From our experiments, it appears that the
model corresponding to maximal degree 3 complex enough to fit the data
well (on the test set, better than other models), but not ``too
complex'' so it does not grossly overfit the training set.

\begin{notes}
Although this experiment suggest 3 as the best choice, in practice a
more careful evaluation would be required. What we really want to
minimize is the generalization error, of which test error is an
estimate; this estimate, though unbiased, is a random variable,
and has nonzero variance. It is conceivable that with a different test
set, another model say of degree 2, would show better test
performance. A typical methodology in practice consists of training
different models on the same training set, choosing one based on a
separate \emph{validation set} (the role played in our experiments by
the ``test set'') and finally evaluating the chosen model on a
separate test set (the last step is missing in our experiments).
\end{notes}
\end{answer}

\problem{15} In real life, the measured values of the target function
(MPG in this case) will be affected by
noise. That is, if the ``true'' value of the function is
$y_{true}=f(\mathbf{x},\mathbf{w}_{true})$, then the value we record (and the one
included in the data set) is $y=y_{true}+\epsilon$, where $\epsilon$
is a random noise with distribution $p(\epsilon)$.
The weights $\hat{\mathbf{w}}$, found with linear
regression from $\mathbf{x}$ and $y$, are an estimate of the true weights $\mathbf{w}_{true}$.  We would
like to find how \emph{biased} is this estimate (more about the bias
of an estimator in the Problem 2). Assuming $\Ep{p(\epsilon)}{\epsilon}
= 0$, what is the bias $\Ep{p(\epsilon)}{\hat{\mathbf{w}}-\mathbf{w}_{true}}$?

\begin{answer}
First, our assumption about the underlying noiseless model here can be
written as $\mathbf{y}_{true}=\mathbf{Xw}_{true}$. Multiplying both
sides by $\left(\mathbf{X}^T\mathbf{X}\right)^{-1}\mathbf{X}^T$, we get
\[
\mathbf{w}_{true}=
\left(\mathbf{X}^T\mathbf{X}\right)^{-1}\mathbf{X}^T\mathbf{y}_{true}.
\]
Now,
\[
\begin{split}
\Ep{p(\mathbf{\epsilon})}{\hat{\mathbf{w}}}=&
\Ep{p(\mathbf{\epsilon})}{\left(\mathbf{X}^T\mathbf{X}\right)^{-1}\mathbf{X}^T\mathbf{y}}\\
=&\Ep{p(\mathbf{\epsilon})}{\left(\mathbf{X}^T\mathbf{X}\right)^{-1}\mathbf{X}^T(\mathbf{y}_{true}+\mathbf{\epsilon})}\\
=&\Ep{p(\mathbf{\epsilon})}{\left(\mathbf{X}^T\mathbf{X}\right)^{-1}\mathbf{X}^T\mathbf{y}_{true}}+\Ep{p(\mathbf{\epsilon})}{\left(\mathbf{X}^T\mathbf{X}\right)^{-1}\mathbf{X}^T\mathbf{\epsilon}}\\
=&\mathbf{w}_{true}+\left(\mathbf{X}^T\mathbf{X}\right)^{-1}\mathbf{X}^T\Ep{p(\mathbf\epsilon)}{\mathbf{\epsilon}}=\mathbf{w}_{true}.
\end{split}
\]
That is, as was stated in the lecture, $\hat{\mathbf{w}}$ is unbiased.
\end{answer}

\end{enumerate}


\pagebreak

\section*{Problem 2: estimation} 
{\bf Reference:} Lecture 1-2; Recitations; Chapter 4.


The distribution of random variables we deal with often belongs to a particular
\emph{parametric family} ,
$P(x\:;\:\mathbf{\theta})$. For instance, the number $h$ of heads in a
series of $n$ coin tosses
has binomial distribution $P(h;p)=\binom{n}{h}p^h(1-p)^{h}$ with
parameter $p$, the probability of a head in a single toss. A Gaussian
univariate distribution has two parameters - the mean $\mu$ and the
variance $\sigma^2$.  
The \emph{likelihood}
$L(x;\mathbf\theta)$ of the data $x$ given the value of the parameter
$\mathbf\theta$ is the value of $P(x|\mathbf\theta)$; the different
terminology and notation here is used to emphasize that we are concerned
with the function of $\theta$, having observed the data, rather than
predicting the value of the data. 




In this problem, we will derive estimators for the parameters of the
\emph{multinomial distribution}. Consider a sequence of $n$ independent
random trials, each having one of $N$ possible outcomes $1,\ldots,N$. In each
trial, the outcome $k$ occurs with probability $\pi_k$. Of course,
since one of $1,\ldots,N$ necessarily occurs, $\sum_k{\pi_k}=1$. 
The random vector $X_1,\ldots,X_N$, where $X_k$ is the count of times
$k$ was the outcome of a trial, has multinomial distribution
\[
P(X_1=x_1,\ldots,X_N=x_N)=\frac{n!}{\prod_{k=1}^Nx_k!}\prod_{j=1}^N\pi_j^{x_j}.
\]
An example of such distribution, for the case of
$N=2$, is the perhaps more familiar binomial distribution - for
instance, the counts of heads and tails in $n$ coin tosses. The
distribution of the counts of 1 through 6 in a series of $n$ rolls of
a die is multinomial as well.
Remember that the values drawn from the multinomial distribution are the
\textbf{counts} of outcomes $1,\ldots,N$ in $n$ trials, and not the
outcomes themselves.

\begin{enumerate}

\item[] {\bf Maximum Likelihood Estimation} We'd like to find the
parameter setting $\hat\theta$ that best ``explains the data''. More
precisely, we need a measure of how well the distribution
$p(\cdot|\theta)$ with different values of $\theta$ fits the observed
data $x$. Since $p(x|\theta)$ is the probability of reproducing the
data $x$ as samples from the distribution, it seems natural to try to
maximize this probability. This gives the \emph{Maximum Likelihood}
estimation criterion for $\theta$:
\[
\hat{\theta}_{ML}=\argmax_\theta L(x; \theta).
\]
Note that likelihood of a set of independently identically distributed
(i.i.d.) random variables is the product of
their likelihoods (by definition of independence):
$L(x_1,\ldots,x_n;\theta)=\prod_{i=1}^nL(x_i;\theta)$. This expression
can be
quite complicated, and make maximization technically challenging.
However, any
monotonically increasing function of the likelihood will have the same
maxima. One such function
is \emph{log-likelihood}  $\log L(x;\theta)$; taking the log turns the
product into a sum, making derivatives significantly simpler. We will
follow this practice of maximizing the log-likelihood instead of likelihood.


\problem{20} Write down the log-likelihood of a set of observed values
from multinomial distribution, and derive the ML estimator for the parameters
$\pi_i,\;i=1,\ldots,N$. 



\begin{advice}
Note that the effective number of parameters (``degrees of freedom'')
of this distribution is $N-1$ rather than $N$ (make sure you
understand why).
\end{advice}
\begin{answer}
We start with writing down the log-likelihood of a set of observations
$X=(x_1,\ldots,x_n)$, given the multinomial parameter vector $\mathbf{\pi}=[\pi_1,\ldots,\pi_N]$:
\[
\log L(X\;;\;\mathbf\pi)=
\log n!-\log \prod_{k=1}^Nx_k!+\sum_{j=1}^Nx_j\log\pi_j
\]

Note that also we refer to $N$ parameters in the above, there are only
$N-1$ degrees of freedom, due to the constraint $\sum_k\pi_k=1$; we
can express one of the parameters in terms of the other,  say $\pi_N=1-\sum_{k=1}^{N-1}\pi_k$. 

Now, we are interested in maximizing log-likelihood as a function of
$\mathbf\pi$, so we can ignore the terms independent of $\mathbf\pi$:
\[
\begin{split}
\hat{\mathbf{\pi}}_{ML}=&\argmax_\mathbf{\pi}\log L(X\;;\;\mathbf\pi)=
\argmax_\mathbf{\pi}\sum_{j=1}^Nx_j\log\pi_j\\
=& \argmax_{\pi_1,\pi_{N-1}}\sum_{j=1}^{N-1}x_j\log\pi_j+x_N\log\left(1-\sum_{i=1}^{N-1}\pi_i\right),
\end{split}
\]
i.e. $\hat{\mathbf{\pi}}_{ML}$ is the solution of the following system of equations:
\[
\left\{\begin{array}{l}
\frac{\partial}{\partial\pi_1}\log L(X\;;\mathbf\pi)=
\frac{x_1}{\pi_1}-\frac{x_N}{1-\sum_{i=1}^{N-1}\pi_i} =
\frac{x_1}{\pi_1}- \frac{x_N}{\pi_N} = 0\\
\vdots\\
\frac{\partial}{\partial\pi_{N-1}}\log L(X\;;\mathbf\pi)=
\frac{x_{N-1}}{\pi_{N-1}}-\frac{x_N}{1-\sum_{i=1}^{N-1}\pi_i} =
\frac{x_{N-1}}{\pi_{N-1}}- \frac{x_N}{\pi_N} = 0\\
\end{array}\right.
\]
From this we can infer that 
\[
\frac{x_i}{\hat\pi_i} = c\;\;\;\mbox{(constant)}
\]
or, alternatively, $\hat\pi_i = x_i/c$, for all $i=1,\ldots,N$. 
Since $\hat\pi_i$ has to sum to one, we see that $c$ has to be
\[
c = \sum_{j=1}^N x_j
\]
This finally gives the result 
\[
\hat\pi_k=x_k/\sum_{j=1}^Nx_j.
\]
That is, the ML estimate of $\pi_k$ is the observed 
frequency $x_k/n$. 

\end{answer}




\problem{5} Let $\hat{\theta}(X_1,\ldots,X_n)$ be an estimator based
on data $X_1,\ldots,X_n$, where the data are viewed as random samples
from $p(\cdot|\theta)$. We can try to assess how well the estimator
can recover the parameters $\theta$. One useful measure is the {\em
bias} of the estimator. This is defined as the expectation
$\E{\hat{\theta}(X_1,\ldots,X_n)-\theta}$, taken with respect to the
true distribution of the data $X_1,\ldots,X_n$, which in this case is
$p(\cdot|\theta)$. The bias measures whether the estimator
systematically deviates from the true parameters $\theta$ that were
used to generate the data. An estimator is called \emph{unbiased} if
its bias is zero. Show that the ML estimator $\hat{{\pi}_k}_{ML}$ is
indeed unbiased in this sense.
\begin{answer}
To show this, we need to find the expected value of the ML estimator,
with respect to the distribution of the observations. Let $o_k^i$,
$i=1,\ldots,n,\;k=1,\ldots,N$  be
random variables
\[
o_k^i=\begin{cases}
1 & \text{if observation $x_i=k$,}\\
0 & \text{otherwise}.
\end{cases}
\]
Note that $\E{o_k^i}=\pi_k$. Now, $\hat{\pi}_k=\sum_{i=1}^no_k^i/n$,
therefore,
\[
\E{\hat{\pi}_k}=\sum_{i=1}^n\E{o_k^i}/n=\sum_{i=1}^n\pi_k/n=\pi_k,
\]
which proves the lack of bias.
\end{answer}




\item[]{\bf Maximum A-Posteriori (MAP) Estimation}
ML estimator treats the parameter $\theta$ as an unknown fixed value, and tries
to find a value which is most likely, based only on the observed
data. A  more sophisticated approach could be to treat  $\theta$ as a
random variable itself, having a distribution $p(\theta)$. We can then seek the value of $\theta$ that is the best
fit both to our prior belief about $\theta$,
expressed by $p(\theta)$, and to the observed data $x$. For
instance, consider a biased coin, and a series of 20 tosses that all
happen to be heads. ML estimator of the parameter $p$ of binomial
distribution  will tell us that $p=1$, i.e. tails will never
occur. This seems unreasonable, however, since we believe that even
the very biased coin would at least rarely produce a tail. 

The \emph{Maximum A-Posteriori} (MAP) estimator allows us to
incorporate such knowledge into the estimate. In MAP, we 
explicitly maximize the posterior probability of $\theta$ given $x$: 
\[
\hat{\theta}_{MAP} = \argmax_{\theta} p(\theta|x).
\]

Now, from the definition of conditional probability, and using Bayes
rule, we have
\[ 
p(\theta|x)=\frac{p(\theta,x)}{p(x)}=\frac{p(x|\theta)p(\theta)}{p(x)}.
\]
The denominator $p(x)$ is independent of $\theta$; it is, in
fact, equal to $\int_\theta p(x,\theta)p(\theta)d\theta$ (why?). Therefore,
it can be ignored in the maximization, and thus the MAP estimator is
\[
\hat{\theta}_{MAP}=\argmax_\theta p(x|\theta)p(\theta).
\]
Note that as opposed to ML, MAP  estimator is somewhat ``subjective'';
its value depends on the particular prior $p(\theta)$ that one
believes is appropriate.
A commonly used prior for the parameters of the  multinomial
distribution is the \emph{Dirichlet} distribution. This is a
distribution over vectors $\mathbf{\pi}$ such that all elements
$\pi_k>0$  and $\sum\pi_k=1$. The Dirichlet distribution can be
written as
\[
D(\mathbf{\pi_1,\ldots,\pi_N})=\frac{1}{Z}\prod_{k=1}^N\pi_k^{\alpha_k-1}
\]
where $\alpha_k$ are the \emph{hyperparameters} (parameters of the
prior), and $Z$ is a normalization constant needed to make this a
valid distribution.


\problem{20} Derive the MAP estimate of the parameters of multinomial
distribution, using the Dirichlet prior.

Hint: write down the term maximized by MAP as the product of
likelihood and prior, take the log, and follow steps similar
to the derivation of the ML estimator. 


\begin{answer}
The MAP estimate, by definition, is
\[
\hat{\mathbf\pi}_{MAP}=\argmax_{\mathbf\pi}\log p(X|\mathbf\pi)p(\mathbf\pi)=
\argmax_{\mathbf\pi}\left(
\sum_{j=1}^Nx_j\log\pi_j+\sum_{j=1}^N(\alpha_j-1)\log\pi_j\right),
\]
and the corresponding optimization constraints are expressed by the system
\[
\frac{x_k+\alpha_k-1}{\pi_k}=\frac{x_N+\alpha_N-1}{\pi_N},\;\;k=1,\ldots,N-1.
\]
Solving this system in a way very similar to our solution for the ML
estimate, we get
\[
\hat\pi_k=\frac{x_k+\alpha_k-1}{\sum_{j=1}^N(x_j+\alpha_j-1)}.
\]
Note that this is equivalent to ML estimate, but taking into account,
in addition to the observations $X$, a hypothetical
additional set of observations, in which the outcome $K$ occurs
$\alpha_K-1$ times.
\end{answer}



\item (\emph{Optional}) We would like the influence of the prior to be in
some sense proportional to our confidence in it. Suggest a way of
achieving this with the Dirichlet prior.

\begin{answer}
We might believe a priori that the parameters $\pi_j$ should be
something like $p_j$, where $p_j\geq 0$, and $\sum_j p_j = 1$. We'd
like to express this prior belief in a Dirichlet distribution and also
specify the strength of our belief. We can do this easily by setting
$\alpha_k = \alpha p_k$, where larger values of $\alpha$ serve to
concentrate the prior distribution around $p_j$ (we would need more
data to overcome such a prior). You can see this from the previous
result by letting $\alpha\rightarrow\infty$. In this case, $\hat\pi_j
= p_j$ as expected. 
\end{answer}

\problem{5} Show that the influence of any fixed prior on the estimate
of multinomial parameters becomes less significant as the amount of
data increases. That is, show that as $K = \sum_k{x_k}$ increases,
$\hat{\pi_k}_{MAP}$ approaches a limit which is independent of the
prior. This is a sensible conclusion since with sufficient data it
shouldn't matter what we thought about the parameter values prior to
seeing the data.

\begin{answer}
First of all, we notice that
\[
\lim_{n\to\infty}\frac{\alpha_k-1}{n+\sum_j\alpha_j-1}=0.
\]
Therefore,
$\lim_{n\to\infty}\hat\pi_k{}_{MAP}=\lim_{n\to\infty}\frac{x_k}{n+\sum_j\alpha_j-1}$.
Now,
\[
\lim_{n\to\infty}\hat\pi_k{}_{ML}/{\hat\pi_k}{}_{MAP}=
\lim_{n\to\infty}\frac{n+\sum_j\alpha_j-1}{n}=1,
\]
that is, as the number of available observations increases, the MAP
estimate approaches the ML estimate, which is of course independent of
the fixed values of $\alpha_j$.  
\end{answer}




\end{enumerate}

\end{document}






