\documentclass[11pt]{article}

\input fullpage.sty

\parindent=0pt
\parskip=0.5em

\let\ul\underline
\let\grad\nabla

\def\Nats{{\cal N}}

\begin{document}

\title{HW 03}
\author{36-711 Fall 2002}
\date{Due October 11, 2002}
\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{enumerate}
\item {\em Paper and pencil, and C.} The program {\tt bisection.c}
  that we discussed in class computes a root of $f(x)=0$ by the
  bisection method, summarized below:
\begin{verbatim}

   given a,b such that f(a)*f(b)<0 :
   loop forever {
     set x = (a+b)/2;
     if f(x) is approx 0, or a-b is approx 0, break out of the loop;
     if f(x)*f(a) > 0, set a = x; else set b = x;
   }
   return x.

\end{verbatim}
\begin{enumerate}
\item Let $x_{n}$ be the $x$ that is computed after $n$ iterations of
  this loop.  Assume $f(x)$ is a continuous function with a unique
  root $x^*$ in the original interval with endpoints $a$ and $b$.
  Show that
\[
        |x_{n+1} - x_{n}| \leq k \cdot |x_n - x_{n-1}|
\]
for some constant $0<k<1$ and all $n$ sufficiently large.  From this,
show that $x_{n}$ converges to $X^*$, and moreover
\[
        |x_{n} - x^*| \leq k^n\cdot (b-a) 
\]
\begin{quote}
{\em Definition:} A sequence $x_n$ is said to converge at rate
(at least) $\alpha$ to its limit $x^*$ if 
\[
        |x_{n} - x^*| \leq k^{n\alpha} \cdot C 
\]
for some constants $0<k<1$, $\alpha>0$, $C>0$.  Note that it is
sufficient to show that
\[
        |x_{n+1} - x^*| \leq k \cdot |x_n - x^*|^{\alpha}
\]
or even
\[
        |x_{n+1} - x^*| \leq k \cdot |x_n - x^*|^{\alpha}
\]
for some constant $k\neq0$ and all $n$ sufficiently large ($k$ is not
necessarily the same in each expression).  If $\alpha=1$ we say the
sequence converges linearly; if $\alpha=2$ we say it converges
quadratically.
\end{quote}
Part (a) shows that the bisection method converges linearly.

\item If you implement {\tt bisection.c} using single precision ({\tt
    float}'s) in C, it should converge in roughly 24 iterations or
    less.  Why?   If you implement it in double precision ({\tt
    double}'s), how many iterations should it converge in?

\item Perhaps an attractive alternative to looking at the midpoint
  $x=(a+b)/2$ each time through the loop, is to draw a straight line
  from $(a,f(a))$ to $(b,f(b))$ and let $x$ be the place where the
  straight line intersects the $x$-axis.  Modify {\tt bisection.c} to
  do this, and compare its performance with the orginal {\tt
    bisection.c} (by running both programs).  Does the new algorithm
  still converge linearly to the root $x^*$?  (You may need to make
  some assumptions about $f'(x)$ to give a precise answer.)
\end{enumerate}



\item {\em C program.} Recall from statistical theory that if
\[
        L(\theta) = L(\theta;y_1,\ldots,y_n) = \prod_{i=1}^n p(y_i|\theta)
\]
is a likelihood for the $\theta$ and the iid data $y_1, \ldots, y_n$,
then a {\em maximum likelihood estimate} (MLE) of $\theta$ is a value
$\hat\theta$ that maximizes $L(\theta)$.

In practice it is usually better to work with the {\em log-likelihood}
\[
        \ell(\theta) = \log L(\theta) = \sum_{i=1}^n \log p(y_i|\theta)
\]
since this is less likely to underflow in floating point
calculations\footnote{It is possible that we would need to add a
  constant to $\ell(\theta)$ to avoid further overflow problems, but
  that is an issue for another day.}, and since $\ell(\theta)$ is maximized
if and only if $L(\theta)$ is maximized.  

We usually find an MLE by setting the derivative of $\ell(\theta)$
equal to zero
\[
        \ell'(\theta) = \sum_{i=1}^n \frac{d}{d\theta}  \log p(y_i|\theta)
        \equiv 0
\]
(this is the {\em likelihood equation}) and solving for $\theta$;
i.e. we want a root of $\ell'(\theta)=0$.

Please do the following:

\begin{enumerate}

\item Modify {\tt root-nr.c} (or write your own program) to use
  Newton-Raphson to calculate $\hat\theta$ for the data
\verb|1.2 2.4 1.3 1.3 0.0 1.0 1.8 0.8 4.6 1.4|
which is to be modeled as iid observations from a Cauchy likelihood, i.e.,
\[
   p(y_i|\theta)
   = \frac{1}{\pi s} \cdot \frac{1}{1+\left(\frac{y_i-\theta}{s}\right)^2}
\]
with $s$ fixed in advance to be either $s=1$ or $s=0.01$.

\item For each of the two values of $s$, use your program to to find
  the MLE, starting from a large variety of initial guesses.  For each
  $s$, describe your results.  What is going on?

\end{enumerate}

[ please turn in your {\tt C} program, results, and conclusions ]


\item {\em Splus or C.}  This method explores the {\em steepest
    ascent/descent} method.
\begin{enumerate}
\item Verify, by differentiating $f(\ul x + t\cdot\ul q)$ with respect to
  $t$, that the derivative in the direction of $\ul q$ is $\ul
  q^T\grad f(\ul x)$, and therefore the direction of maximum increase
  of $f(\ul x)$ is $\ul q = \grad f(\ul x)$, the gradient of $f(\ul
  x)$. [and the direction of maximum {\em decrease} is $\ul q = -\grad
  f(\ul x)$].

\item Write a routine in Splus that implements the ``method of
  steepest descent'' (Monahan, p. 188), and applies it to the function
  in Example 8.4 (same page of Monahan).  The basic algorithm is
\begin{verbatim}
given subroutines that calculate f(x) and df(x)=gradient of f, 
and given starting value x0:

 loop forever {

  let x0=xnew
  let u minimize f(x0 + u*grad f(x0)) 
  let xnew=x0 + u*grad f(x0)
  if relative change from x0 to xnew is small, or xnew is near zero
      break out of the loop
}
\end{verbatim}
Notes: 
\begin{itemize}
\item In Splus you can use {\tt nlmin} to do the minimization, for
    example, or code one of the univariate search methods by hand.  In
    C you have these minimization routines already, but you would have
    to code evaluation of the vector functions by hand.
\item If it seems too painful to code a search method, one could
  simply start with a stepsize of $u=1$, and repeatedly half the
  stepsize until $f(x_0 + u\cdot\grad f(x_0)) < f(x_0)$, for example.
  If you try this, how does this compare with using the optimal $u$?
\end{itemize}
Verify (and plot) the zigzag behavior illustrated in the figure on
p. 189 of Monahan.

\item Refer to the modification of steepest descent by ``taking a step
  in the Newton direction'' at the bottom of p. 188.  Do Ex. 8.19, p.
  197 of Monahan, to obtain the $u$ that minimizes the quadratic
  approximation described on p. 197.

\item  Use your result in part (b) to revise your routine from part
  (a).  Check (and plot) the behavior of the new algorithm.  Does it
  still converge to the right answer?  Does it zigzag?  Etc. \\ 
Notes:
\begin{itemize}
\item Monahan states on p. 188 that the relevant step comes from
  equation (8.5.2) but in fact it comes from equation (8.5.4), so you
  will need a subroutine that can compute the Hessian as well.
\item Again, an alternative to calculating $u$ according to the recipe
  of Ex. 8.19 would be to start with $u=1$ and repeately half it until 
  $f(x_0 + uH^{-1}(\ul x_0)\cdot\grad f(x_0)) < f(x_0)$.
\end{itemize}

\end{enumerate}

\end{enumerate}

\end{document}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Great exercise, requires update of ``calling C from S'' though...

\item {\em Paper and pencil, Splus and C.} Following are three 
  Splus subroutines {\tt nll1.s}, {\tt nll2.s} and {\tt nll3.s}, that
  each compute the log-likelihood for iid N(0,1) observations [these
  are available in the class data directory].

\begin{verbatim}
-----------------------------------------------------------------
nll1 _ function(x,theta) {
#
# nll1
# 
# compute the normal(theta,1) loglikelihood for the data in the vector
# x, at each of the values in the vector theta
#
        n _ length(x)
        m _ length(theta)
        y _ rep(0,m)

        for (i in 1:m) {
           y[i] _ -n/2 * log(2*pi)
           for (j in 1:n) 
                y[i] _ y[i] - (x[j] - theta[i])^2/2
           }

        return(y)

}

---------------------------------------------------------------

nll2 _ function(x,theta) {
#
# nll2
# 
# compute the normal(theta,1) loglikelihood for the data in the vector
# x, at each of the values in the vector theta
#
        n _ length(x)
        m _ length(theta)
        y _ rep(0, m)
        for(i in 1:m) {
                y[i] _  - n/2 * log(2 * pi) - sum((x - theta[i])^2)/2
        }
        return(y)

}

---------------------------------------------------------------

nll3 _ function(x,theta) {
#
# nll3
# 
# compute the normal(theta,1) loglikelihood for the data in the vector
# x, at each of the values in the vector theta
#
#
# NOTE: The following was how it was originally coded, but this is 
# not as fast as nll2, and it's not much easier to read!
#
#       n _ length(x)
#       m _ length(theta)
#
#       delta _  outer(x,theta,"-")
#       SS _ apply(delta^2/2,2,sum)
#       y _ -n/2 * log(2*pi) - SS/2
#
#       return(y)
#

  - length(x)/2 * log(2 * pi) - apply(outer(x, theta, "-")^2/2, 2, sum)

}

-------------------------------------------------------------
\end{verbatim}

For example, you could plot a normal log-likelihood for 100 N(0,1)
observations like this:
\begin{verbatim}
> x _ rnorm(100)                # or some real data if you have it
> th _ seq(-3,3, length=100)
> X11()                         # or motif() if you prefer
> plot(th,nll3(x,th),type="l")
\end{verbatim}
\begin{enumerate}

\item The guts of {\tt nll3.s} are the lines
\begin{verbatim}
        delta <- outer(x,theta,"-")
        SS <- apply(delta^2,2,sum)
        y <- -n/2 * log(2*pi) - SS/2
\end{verbatim}
  Explain in detail what is going on in these lines.  Use a small
  example, e.g. use explicit vectors of length 2 or 3 for both {\tt x}
  and {\tt theta}.

\item Write a new normal log-likelihood function that looks something
like this:
\begin{verbatim}
nll4 <- function(x,theta) {
#
# nll4
# 
# compute the normal(theta,1) loglikelihood for the data in the vector
# x, at each of the values in the vector theta.  Do all the
# computations in C.
#
        n <- length(x)
        m <- length(theta)
        y <- rep(0,m)

        z <- .C("mynll",as.integer(n),as.integer(m),
                 as.double(x),as.double(theta),ans = as.double(y))

        return(z$ans)
     
\end{verbatim}
Of course you will also have to write a {\tt C} subroutine {\tt mynll}
to do the computation, and load it into Splus using {\tt dyn.load2()}.

Which of the four routines  {\tt nll1}, {\tt nll2}, {\tt nll3} and {\tt
nll4}  seems to be fastest?

\item To get a better answer to the last question in part (b), write
a function like the following in Splus
\begin{verbatim}
nlltime <- 
function(nll, x, th, n = 100)
{
# nlltime
        print(unix("date"))
        for(i in 1:n)
                nll(x, th)
        print(unix("date"))
        invisible()
}
\end{verbatim}
Now the command {\tt nlltime(nll2,x,theta)} will show approximately
how long it takes to run {\tt nll2} 100 times; similarly for {\tt
nlltime(nll3,x,theta)}, etc.

Experiment with different lengths of {\tt x} and {\tt theta} to find
out which of the four nll routines is fastest under which conditions.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
