\documentstyle[11pt]{article}
\setlength{\parindent}{0cm}
\setlength{\parskip}{0.2cm}
\textwidth 6in
\textheight 8.5in
\topmargin -0.375in
\oddsidemargin 0.375in
\newcommand{\fol}{\mbox{$\prec \prec$}}
\newcommand{\iuc}{\mbox{${\cal I}^c$}}
\newcommand{\ilc}{\mbox{${\cal I}_c$}}
%\newcommand{\qued}{\rule{2mm}{3.5mm}}
\input{extras}
\begin{document}
\mbox{} \hfill {\sf Draft: \today}\\[0.2cm]
% Draft 0. (???)
\begin{centre}
{\Large \bf An Algorithm for Unimodal Isotonic Regression, with
Application to Locating a Maximum}\\[0.2cm]
by\\[0.2cm]
R. A. Mureika and T.R. Turner\\[0.2cm]
{\em Department of Mathematics and Statistics,\\[0.2cm]
University of New Brunswick}\\[0.2cm]
and\\[0.2cm]
P. C. Wollan\\[0.2cm]
{\em Division of Biostatistics,\\[0.2cm]
Mayo Clinic}\\[0.5cm]

{\bf Abstract}\\
\end{centre}

A new algorithm is proposed for applying isotonic regression to data
having an underlying unimodal structure.  This algorithm is actually
a special case of a more general procedure which permits a ``divide
and conquer'' approach to a broad class of isotonic regression
problems.  It is noted that repeated application of the algorithm
permits the estimation of the location of the maximum of a data set
having underlying unimodal means, and that this location estimate is
consistent.  The performance of the resulting procedure for locating
a maximum is assessed through a simulation study.

{\small Key words: {\em unimodal ordering, level sets,
level values, divide and conquer approach, tree-like
orderings}.}\\[0.2cm]

\section{Introduction}
\label{intro}
Algorithms for implementing isotonic regression under orderings other
than the simple linear order are difficult to construct.
The best known of such algorithms is the Maximum Lower Sets algorithm
(Robertson, Wright, and Dykstra, 1988, p.24);
this is complicated and hard to program.  It is
also reputed to run rather slowly, and indeed the number of operations
required grows exponentially in certain cases.

The motivation for developing an improved algorithm for performing
such regressions came in part from a data set being studied by
members of the Faculty of Forestry at the University of New Brunswick.
These data consisted of observations which had been made of the
``vigour'' of growth of five stands of black spruce.  The stands each
had different initial tree densities.  It was expected that vigour
would initially increase (as the trees increased in size) and then
level off and start to decrease as the growing trees encroached upon
each others' space and competed more strongly for resources such as
moisture, nutrients, and light.  It was further expected that the
position of the mode of the vigour observations would depend upon the
initial densities.

Plots of the data did not make it completely clear as to where the
leveling-off point or mode occurred; the Forestry researchers
requested a procedure for determining the location of this mode.  A
procedure which comes immediately to mind is to fit unimodal isotonic
regressions with mode at each of the possible locations in turn.  The
location yielding minimal error sum of squares is then chosen as the
location of the mode.  It is thus desirable to be able to perform a
large number of unimodal isotonic regressions quickly and
efficiently.

Formally the unimodal isotonic regression problem may be stated as
follows:  Suppose that $Y_{ij}$, $i=1, \ldots, p$, $j = 1, \ldots,
n_i$, are independent random variables such that $Y_{ij} = \mu_i +
E_{ij}$ for all $i$ and $j$, where the $E_{ij}$ have mean 0 and
variance $\sigma^2$.  Further suppose that the $\mu_i$ have a
{\em unimodal ordering}, i.e. that

\begin{equation}
\label{unimod1}
\mu_1 \leq \mu_2 \leq \ldots \leq \mu_{k_0} \geq \mu_{k_0 + 1}
\geq \ldots \geq \mu_p \;.
\end{equation}

It is desired to estimate the values of $\mu_1, \ldots,
\mu_p$.  The (weighted) least squares estimates of the
$\mu_i$ are given by minimizing
\[
SS = \sum_i \sum_j (Y_{ij} - \hat{\mu}_i)^2 w_i
\]
subject to the constraint (\ref{unimod1}), where $w_1,
\ldots, w_p$ are a (given) set of positive weights.  This
is of course a particular case of isotonic regression under
a partial ordering.  It may initially be subdivided into
three sub-problems involving only {\em linear} orderings:
(a) estimating $\mu_1, \ldots , \mu_{k-1}$; (b) estimating
$\mu_k$; and (c) estimating $\mu_{k+1}, \ldots , \mu_n$.
Sub-problem (b) is totally trivial and sub-problems (a) and
(c) can be solved by standard and well-known techniques.
The question is how to combine the solutions of the three
subproblems appropriately.

The answer is essentially to ``interleaf'' the estimates
resulting from solving sub-problems (a) and (c) in {\em
numerical} order, tack on $\hat{\mu}_k = \bar{Y}_{k.}$ at the upper
end, solve the corresponding isotonic regression with
respect to the resulting linear ordering, and then put the
estimates back in their original order.

In the next section some notation and terminology are established,
and the main result, which validates the proposed algorithm, is stated.
In the following section, \ref{generalize}, it is indicated how this
result may be applied to other ``tree-like'' orderings as well as to
unimodal ones.  In section \ref{locmax} the problem of estimating the
location of a maximum is discussed in more detail.  The procedure of
performing unimodal isotonic regressions for all possible mode
locations is shown to yield a consistent estimate of the true mode
location.  Simulation studies assessing the procedure are then
discussed.  Section \ref{concl} consists of concluding remarks.
Some technical lemmas, upon which the proof of the main result depends,
their proofs, and the proof of the main result itself, are
given in an appendix.

\section{Notation and Terminology, and the Main Result}
\label{mainres}
Let $S$ be a set and let $\prec$ be a partial order on
$S$.  Recall that an isotonic function (with respect to the
partial order $\prec$) is a (real-valued) function $f$ such
that $x \prec y$ implies $f(x) \leq f(y)$.  If $g$ is an
arbitrary function on $S$, and $w$ is a non-negative
(weight) function on $S$, then the {\em isotonic
regression} of $g$, with respect to $\prec$ and $w$,
(denoted $g_*$) is that value of $\hat{g}$ which minimizes
\begin{displaymath} \sum_{s \in S} [g(s) -
\hat{g}(s)]^2w(s) \end{displaymath} over all {\em isotonic}
functions $\hat{g}$.

Let $S_1$ and $S_2$ be two subsets of $S$.  We say that $S_2$ {\em
follows} $S_1$, (in symbols $S_1 \fol S_2$) if $x \prec y$ for every
$x$ in $S_1$ and every $y$ in $S_2$.

Now let $S = S_1 \cup S_2$ where $S_1$ and $S_2$ are
disjoint, and $S_1 \fol S_2$.  Let $g_j$ be the restriction
of $g$ to $S_j$, and suppose that $g_{j*}$ is the isotonic
regression of $g_j$, $j = 1,2$.  The weight functions used
to form $g_{j*}$ above are of course the restrictions of
the overall weight function $w$.

An elementary but important fact about isotonic regression
is that $g_*$ takes the form \begin{displaymath} g_*(s) =
c_i \mbox{ on } L_i, \; i = 1, \ldots, r \end{displaymath}
where $L_1, \ldots, L_r$ form a disjoint and exhaustive
collection of subsets of $S$, and $c_1 < c_2 < \ldots <
c_r$.  Moreover $c_i$ is the weighted mean over $L_i$ of
the values of $g(s)$; i.e.  \begin{displaymath} c_i =
\frac{\sum_{s \in L_i} w(s)g(s)}{\sum_{s \in L_i}
w(s)}\;\;.  \end{displaymath} (See
Robertson, Wright, and Dykstra, 1988, p. 18, Theorem 1.3.5.)
We call the sets $L_i$ the {\em level} sets
and the values $c_i$ the {\em level} values of the isotonic
regression.

Let the level sets and
level values for $g_{1*}$ be $L_1, \ldots, L_{r_1}$ and
$c_1 <  \ldots <  c_{r_1}$, and those for $g_{2*}$ be
$L_{r_1+1}, \ldots, L_r$ and $c_{r_1+1} <  \ldots <  c_r$.
(Note:  It may well occur that $c_{r_1} \geq c_{r_1+1}$.)
Define a function $f$ on $\{1, \ldots, r\}$, by $f(t) =
c_t$ for $t = 1, \ldots r$, and a weight function $u$ by
\begin{displaymath} u(t) = \sum_{x \in L_t} w(x) \;\;.
\end{displaymath}

We are now ready to state:

{\bf Theorem 1:} Let $S_1$, $S_2$, $f$, and $u$ be as in
the foregoing discussion.  Let $f_*$ be the isotonic
regression of $f$ with respect to the usual order on $\{1,
\ldots, r\}$ and the weight function $u$.  Then the
isotonic regression of $g$ with respect to $\prec$ and $w$
is given by \begin{displaymath} g_*(s) = f_*(t) \mbox{ for
} s \in L_t \;\;.  \end{displaymath}

Remark:  In the case of the unimodal structure, the r\^{o}le of $S_2$
is played by the singleton $\{k\}$, and that of $S_1$ by $\{1,
\ldots, k-1, k+1, \ldots, n \}$.  For the partial order determining
the structure of interest, the two sets $\{1, \ldots, k-1 \}$ and
$\{k+1, \ldots, n \}$ are unrelated.  It is therefore easy to see,
(and well-known; see, e.g., Robertson, Wright, and Dykstra, 1988,
p.57), that an isotonic regression on their union is simply the
amalgamation of separate isotonic regressions on each component.

\section{Generalization}
\label{generalize}
Theorem 1 permits a ``divide and conquer'' approach to
isotonic regression problems.  If the set $S$ decomposes
into two disjoint subsets $S_1 \fol S_2$, and if the
restricted isotonic regression problem can be solved for
each $S_i$, then the overall problem can be solved.
Applying this notion inductively we immediately obtain the
generalization of Theorem 1 to sequences of more than two
sets:

{\bf Corollary 1:} If $S = \cup_{k=1}^n S_k$ and $S_1 \fol S_2 \fol
\ldots \fol S_n$, and the $S_i$ are pairwise disjoint, then the
isotonic regression on $S$ can be calculated by first calculating the
isotonic regression separately on each of $S_1, S_2, \ldots , S_n$,
placing the resulting level values in order and carrying out the
weighted isotonic regression of these level values with respect to
simple linear order.  The weights in question are \[ \sum_{s \in L_i}
w(s) \] where $L_i$ is the level set corresponding to a level value
$c_i$.

The ``divide and conquer'' approach may also be applied to
more complicated arrangements amongst the set $S_i$, called
``tree-like'' structures.\\
{\bf Definition:} Suppose that $S = \bigcup_{k=1}^n S_k.$
We say that $S_1, S_2, \ldots, S_n$ have a tree-like
structure with respect to $\fol$ if \begin{list}{}{}
\item[(i)] the $S_i$ are pairwise disjoint, \item[(ii)] no
element of $S_i$ is comparable with any element of $S_j$
for $i \neq j$ unless $S_i \fol S_j$ or $S_j \fol S_i$, and
\item[(iii)] for any $S_j$ there is at most one $S_i$ with
the property that (1) $S_j \fol S_i$ and (2) there is {\em
no} $S_k$ for which $S_j \fol S_k \fol S_i$.  \end{list}

Now suppose that $S_1, S_2, \ldots, S_n$ have a tree-like
structure, and that the restricted isotonic regression
problem can be solved for each $S_i$.  Then the isotonic
regression on $S$ can be carried out thusly:  Group those
$S_j$ which follow {\em no sets} into collections having a
common set by which they are all followed.  Denote any one
such collection by $S_{i_1}, S_{i_2}, \ldots, S_{i_r}$ with
predecessor $S_i$.  Let $S_1' = S_{i_1} \cup \ldots \cup
S_{i_r}$, and let $S_2' = S_i$.  Since there is no order
relation amongst the elements of the sets $S_{i_j}$, the
isotonic regression on their union is the simple
amalgamation of the individual isotonic regressions.  Thus
the restricted isotonic problem can be solved on $S_1'$; by
assumption it can be solved on $S_2'$.  Therefore, by
Theorem 1 it can be solved on $S_1' \cup S_2'$.  We now
remove $S_i$ and $S_{i_1}, \ldots, S_{i_r}$ from the
collection $S_1, \ldots, S_n$ and replace them by their
union $=S_1' \cup S_2'$.  This yields a new, smaller,
tree-like structure on each component of which the
restricted isotonic problem can be solved.  We now repeat
the procedure until the complete problem is solved, i.e.
until the resulting tree-like structure consists of the
single set $S$.

\section{Estimating the Location of a Maximum}
\label{locmax}
Let $Y_{ij}$ and $w_i$, $i=1, \ldots, p$, $j = 1, \ldots, n_i$, be as
described in section \ref{intro}.  Suppose that the value of $k_0$ is
unknown and one wishes to estimate it in some rational manner.  The
(weighted) least squares estimate of $k_0$ may be determined by
assuming that $k_0 = k$ for each $k = 1, \ldots , p$ and finding the
(weighted) least squares estimates of the $\mu_i$, say
$\hat{\mu}_i(k)$ under this assumption.  Let $SS(k)$ be the
corresponding error sum of squares, i.e.
\[
SS(k) = \sum_i \sum_j (Y_{ij} - \hat{\mu}_i(k))^2 w_i
\]
The estimated value of $k_0$ is then that value of k which minimizes
$SS(k)$.

If we assume that the mode is a strict one, i.e. that
\begin{equation}
\label{unimod2}
\mu_1 \leq \mu_2 \leq \ldots \leq \mu_{k_0 - 1} < \mu_{k_0} > \mu_{k_0 + 1}
\geq \ldots \geq \mu_p \;,
\end{equation}
then it is not hard to demonstrate that this procedure yields
a consistent estimate of $k_0$:

Note that
\[
SS(k) = \sum_i \sum_j (Y_{ij} - \hat{\mu}_i(k))^2 w_i =
        \sum_i \sum_j (Y_{ij} - \bar{Y}_{i.})^2 w_i +
        \sum_i (\bar{Y}_{i.} - \hat{\mu}_i(k))^2 n_i w_i\;.
\]

The vector $[ \hat{\mu}_i(k) ]_{i=1}^p$ is the weighted least squares projection
% \hat{\mu}_i(k)
(with weights $n_i w_i$) of the data onto the cone determined by the hypothesis
\begin{equation}
\label{cone}
\mu_1 \leq \mu_2 \leq \ldots \leq \mu_{k} \geq \ldots \geq \mu_p \;.
\end{equation}
As each $n_i \rightarrow \infty$, $\bar{Y}_{i.} \rightarrow \mu_i$,
and $\hat{\mu}_i(k)$ converges to  $\tilde{\mu}_i$ where
$[\tilde{\mu}_i]_{i=1}^p$ is the isotonization
of $[\mu_i]_{i=1}^p$, that is
% $[\mu_i]$
the weighted least squares projection of $[\mu_i]_{i=1}^p$
onto the cone (\ref{cone}),
with weights
\[
u_i = \lim_{\mbox{ all } n_i \rightarrow \infty} \frac{n_i w_i}{\sum_j n_j w_j}
\]

If $k \neq k_0$, then $\mu_i \neq \tilde{\mu}_i$ for at
least one $i$, and so $SS(k) \geq (\bar{Y}_{i.} -
\hat{\mu}_i(k))^2 n_i w_i \sim$ \linebreak $(\mu_i -
\tilde{\mu}_i)^2 n_i w_i \rightarrow \infty$.  If $k$ does
equal $k_0$, then the distribution of $SS(k)$ is
asymptotically $\sigma^2 \bar{\chi}^2$ (Shapiro, 1985) and so
$P[SS(k_0) \geq SS(k)] \rightarrow 0$ for $k \neq k_0$.
Hence if $\hat{k}$ is the value that minimizes $SS(k)$,
then $P[\hat{k} \neq k_0] \rightarrow 0$.

The isotonic approach thus gives a consistent estimator of
the location of the mode of a unimodal set of means.  The
algorithm described in sections \ref{intro} and
\ref{mainres} provides an effective procedure for
calculating this estimator.

The sequence of means (\ref{unimod2}) could be thought of
as values along a regression curve.  The classical approach
to locating a maximum, dating back at least to Hotelling (1941)
is to fit a parametric model, often a
low-degree polynomial, and to find the maximum value of the
fitted curve.  More recently, nonparametric regression
methods have been applied.  Silverman (1985)
proposed the use of smoothing splines.  Kernel estimators
have been used more for the related problem of estimating
the mode of a probability density (Parzen, 1962,
and, for example, Romano, 1988) but have also
been used for locating maxima (M\"{u}ller, 1989)

Both the parametric and nonparametric approaches have
serious drawbacks:  the parametric approach requires
careful choice of the form of the function, while
nonparametric regression requires choice of a smoothing
parameter.  The method proposed above requires neither.  In
some ways, it is more similar to ranking and selection
methods (see, for example Gibbons, 1977)
than to regression, in that the maximum is chosen to be one
of the observed values.  However, the isotonic approach
exploits the unimodal structure of the means to yield
additional power.

To investigate the efficacy of the isotonic approach, and
to compare this approach with three other possible methods,
a set of simulations was undertaken.  The three other
methods considered were  \begin{list}{}{} \item[(i)]
Fitting (by least squares) a (concave downward) quadratic
and choosing the largest fitted value.  \item[(ii)] Fitting
(by least squares) a (concave downward) ``two-stick''
model, with the bend constrained to come in turn at each of
the possible modal values.  The modal value giving the
minimal sum of squares is then chosen as the maximum
point.  \item[(iii)] Simply picking the largest data
value.  \end{list}

Data to test all four methods were generated from three underlying
models:  \begin{list}{}{} \item[(a)] A (concave downward) quadratic,
with maximum at $x = k$ and height 1 above the $x$-axis, and
$x$-intercepts a distance $n-1$ apart (data being generated at points
$1, 2, \ldots, n$).  \item[(b)] A (concave downward) ``two-stick''
model, with maximum at $x = k$ and height 1 above the $x$-axis and
intercepts at $0$ and $n$.  \item[(c)] A delta function equal to 1 at
$x = k$ and equal to 0 elsewhere on $[1,n]$.  \end{list}

To each model was added pseudo-random independent Gaussian noise,
with zero mean and constant standard deviation $\sigma$.  Six
different choices of the triple $(n,k,\sigma)$ were investigated,
giving a total of 18 simulations, which are tabulated as follows in
Table 1.

\begin{table}[htb]
\begin{centre} {\bf Table 1:  Structure
of Simulations}\\[0.5cm]

\begin{tabular} {| r | r | r | r | r |} \hline
sim. no. & model & $n$ & $k$ & $\sigma$ \\ \hline
1 & 2-stick     & 19          & 10          & 0.5 \\ \hline
2 & quad.       & 19          & 10          & 0.5 \\ \hline
3 & delta fn.   & 19          & 10          & 0.5 \\ \hline
4 & 2-stick     & 19          & 5           & 0.5 \\ \hline
5 & quad.       & 19          & 5           & 0.5 \\ \hline
6 & delta fn.   & 19          & 5           & 0.5 \\ \hline
7 & 2-stick     & 50          & 12          & 0.5 \\ \hline
8 & quad.       & 50          & 12          & 0.5 \\ \hline
9 & delta fn.   & 50          & 12          & 0.5 \\ \hline
10 & 2-stick     & 19          & 10          & 1   \\ \hline
11 & quad.       & 19          & 10          & 1   \\ \hline
12 & delta fn.   & 19          & 10          & 1   \\ \hline
13 & 2-stick     & 19          & 5           & 1   \\ \hline
14 & quad.       & 19          & 5           & 1   \\ \hline
15 & delta fn.   & 19          & 5           & 1   \\ \hline
16 & 2-stick     & 50          & 12          & 1   \\ \hline
17 & quad.       & 50          & 12          & 1   \\ \hline
18 & delta fn.   & 50          & 12          & 1   \\ \hline
\end{tabular} 
\end{centre}
\end{table}

Five hundred replications were generated from each model,
with each of the foregoing parameter combinations.  For
each replication, the location of the mode was estimated
using each of the four possible techniques.  To summarize
the results, the percentage of correct answers and the
mean-squared error of the estimates, over the five hundred
replications, was calculated for each estimation method.
The results are given in Tables 2 and 3.

\begin{table}[htb]
\begin{centre} 

{\bf Table 2:  Percentage Correct}\\[0.5cm]

\begin{tabular} {| r | r | r | r | r |} \hline
         & \multicolumn{4}{c|}{ estimation method } \\ \hline
sim. no. & ``pick max.'' & quad. fit & ``2-stick'' fit & isotonic \\ \hline
1 & (4) 18.6   & (1) 42.2   & (3) 19.8   & (2) 21.0   \\ \hline
2 & (3) 11.8   & (1) 42.2   & (4) 9.4   & (2) 13.2   \\ \hline
3 & (2) 57.2   & (4) 10.4   & (3) 21.2   & (1) 58.4   \\ \hline
4 & (3) 19.0   & (4) 2.8   & (1) 25.8   & (2) 19.8   \\ \hline
5 & (3) 14.4   & (2) 19.8   &  (4) 7.2   & (2) 14.6   \\ \hline
6 & (1) 54.4   &  (4) 3.2   & (3) 23.0   & (1) 54.4   \\ \hline
7 &  (3) 9.8   &  (4) 0.4   & (1) 14.8   & (2) 11.2   \\ \hline
8 &  (3) 4.4   & (1) 12.2   &  (4) 2.0   &  (2) 5.2   \\ \hline
9 & (1) 41.0   &  (4) 0.2   &  (3) 9.6   & (2) 40.2   \\ \hline
10 & (3) 12.6   & (1) 21.0   &  (4) 9.6   & (2) 13.0   \\ \hline
11 & (3) 10.6   & (1) 19.4   &  (4) 6.0   & (2) 11.0   \\ \hline
12 & (1) 26.6   & (4) 10.4   & (3) 11.0   & (2) 26.0   \\ \hline
13 & (2) 11.6   &  (4) 3.4   & (3) 11.2   & (1) 12.0   \\ \hline
14 &  (4) 7.4   &  (2) 8.0   &  (1) 8.2   &  (2) 8.0   \\ \hline
15 & (2) 23.2   &  (4) 1.0   & (3) 12.8   & (1) 24.2   \\ \hline
16 &  (1) 6.4   &  (4) 1.4   &  (4) 6.2   &  (1) 6.4   \\ \hline
17 &  (3) 3.2   &  (2) 3.6   &  (4) 3.0   &  (1) 3.8   \\ \hline
18 & (1) 11.2   &  (4) 1.8   &  (3) 5.0   & (2) 10.8   \\ \hline
\end{tabular} 
\end{centre}

\mbox{}\\ {\footnotesize Remark:  Numbers in parentheses
indicate the rank of the method among the four methods,
according to the ``percentage correct'' criterion.}

\end{table}

\begin{table}[htb]
\begin{centre}
{\bf Table 3:  Mean Square Error}\\[0.5cm]

\begin{tabular} {| r | r | r | r | r | r | r | r |} \hline
         & \multicolumn{4}{c|}{ estimation method } \\ \hline
sim. no. & ``pick max.'' & quad. fit & ``2-stick'' fit & isotonic \\ \hline
1 &  (3) 8.4   &  (1) 1.8   & (4) 10.5   &  (2) 5.8   \\ \hline
2 &  (3) 9.9   &  (1) 2.3   & (4) 16.2   &  (2) 9.1   \\ \hline
3 & (2) 13.7   & (3) 14.4   & (4) 34.1   & (1) 11.9   \\ \hline
4 & (4) 16.6   & (3) 13.8   & (1) 10.9   & (2) 12.6   \\ \hline
5 &  (3) 9.0   &  (1) 4.5   & (4) 13.9   &  (2) 8.4   \\ \hline
6 & (2) 25.9   & (3) 31.1   & (4) 31.9   & (1) 23.9   \\ \hline
7 &  (4) 82.3  &  (3) 78.9  &  (1) 30.2  &  (2) 52.8  \\ \hline
7 &  (3) 50.6  &  (1) 16.9  &  (4) 98.3  &  (2) 39.7  \\ \hline
9 & (2) 246.5  & (3) 275.9  & (4) 346.0  & (1) 231.6  \\ \hline
10 & (3) 16.0   &  (1) 7.6   & (4) 25.2   & (2) 13.7   \\ \hline
11 & (3) 16.9   &  (1) 8.9   & (4) 28.0   & (2) 15.5   \\ \hline
12 & (3) 24.4   & (1) 18.6   & (4) 40.4   & (2) 23.0   \\ \hline
13 & (3) 28.4   & (1) 22.7   & (4) 28.7   & (2) 26.0   \\ \hline
14 & (3) 13.2   &  (1) 7.8   & (4) 19.3   & (2) 11.4   \\ \hline
15 & (4) 41.4   & (3) 38.0   & (1) 37.8   & (1) 37.8   \\ \hline
16 & (4) 186.7  & (1) 108.9  & (3) 124.4  & (2) 117.6  \\ \hline
17 &  (3) 96.0  &  (1) 35.2  & (4) 113.3  &  (2) 69.4  \\ \hline
18 & (3) 339.8  & (2) 324.7  & (4) 362.0  & (1) 299.3  \\ \hline
\end{tabular} 
\end{centre}

\mbox{}\\ {\footnotesize Remark:  Numbers in parentheses
indicate the rank of the method among the four methods,
according to the ``mean square error'' criterion.}

\end{table}

Histograms of the estimation errors were plotted for each
experiment.  An illustrative example of these (for
simulation 13) is given in Figure 1.

\begin{centre} (Figure 1 about here.) \end{centre}

The bias of the quadratic method, toward the centre of the
set of means when the location of the maximum is
off-centre, is readily apparent in this illustration.

The isotonic technique performed consistently well.  It was
either best or second best of all four methods in all
cases, according to both of the above criteria.  It was
outperformed by the quadratic estimation method when the
underlying model was quadratic and the noise level was
low.  With the higher noise level, the isotonic method beat
or tied quadratic estimation in two of the three cases.
The quadratic estimation procedure also did well when the
underlying model was the ``two-stick'' one, particularly
when the maximum was at the centre of the sequence of
means.  When the underlying model was the delta function,
the ``pick the max'' procedure was always the nearest
competitor for the isotonic method.  In fact, the isotonic
method and the ``pick the max'' method tended to give the
same answer a large percentage of the time:  Between 72.0\%
(simulation 18) and 94.2\% (simulation 5).  The
``two-stick'' estimation procedure performed well
essentially only when the underlying model was the
``two-stick'' one.  For simulation number 14, with an
underlying quadratic model, it came first according to the
percentage correct criterion, but last according to the
mean square error criterion.

\section{Concluding Remarks}
\label{concl}
The proposed algorithm for performing unimodal isotonic
regression is simple to understand and to program, and
performs satisfactorily in practice.  It makes feasible the
procedure for estimating the location of a maximum by
performing unimodal isotonic regressions for all possible
locations.  Our simulation studies indicate that this
technique for estimating the position of a maximum is a
sound procedure, over-all.  It works effectively
irrespective of the nature of the underlying unimodal
structure of the means.

The results established in section 3 have of course much
more general application than to the problem of estimating
the location of a maximum.

The software to implement all of the relevant computations,
including the simulations, was written in the S programming
language.  This software is available from the second
author upon request.

{\bf Acknowledgement:}  The authors would like to thank
Kirk Schmidt, a graduate student in the Department of
Forest Engineering, U.N.B., and his advisor Professor Ted
Needham, for drawing the problem on tree growth vigour to
their attention.

\section{Appendix: Proof of the Main Result}

{\bf Definition:} For any constant $c$ we define
\begin{displaymath} \iuc = \{g | g \mbox{ is isotonic and }
g(s) \leq c \mbox{ for all } s \in S \} \end{displaymath}
and \begin{displaymath} \ilc = \{g | g \mbox{ is isotonic
and } g(s) \geq c \mbox{ for all } s \in S \}
\end{displaymath}

Let $g_*(s)$ be the isotonic regression of $g$ and define
\begin{displaymath} g_{cu}(s) = \left \{ \begin{array}{cl}
g_*(s) & \mbox{ if } g_*(s) \leq c\\ c & \mbox{ if } g_*(s)
> c \;\;.\end{array} \right.  \end{displaymath}

{\bf Lemma 1:} The function $g_{cu}$ uniquely minimizes
\begin{equation} \sum_{s \in S} [g(s) - \hat{g}(s)]^2 w(s)
\label{eq:trunciso} \end{equation} subject to $\hat{g} \in
\iuc$.

{\bf Proof:} For any $\hat{g}$ in $\iuc$, \begin{eqnarray*}
\sum_{s \in S} [g(s) - g_{cu}(s)][g_{cu}(s) - \hat{g}(s)
]w(s) & = & \sum_{s \in S} [g(s) - g_*(s)][g_{cu}(s) -
g_*(s)]w(s)\\
  & & + \sum_{s \in S} [g_*(s) - g_{cu}(s)][g_{cu}(s) -
					     \hat{g}(s)]w(s)\\
  & & + \sum_{s \in S} [g(s) - g_*(s)][g_*(s) -
  \hat{g}(s)]w(s)\\ & = & \Sigma_1 + \Sigma_2 + \Sigma_3
\end{eqnarray*} Now $ \Sigma_1 = 0 $ by Theorem 1.3.6 of
Robertson, Wright, and Dykstra, (1988, p.21),
since $g_{cu}(s) - g_*(s)$ is a function
of $g_*(s)$; $ \Sigma_3 \geq 0 $ since $g_*$ is the
isotonic regression of $g$ (applying Theorem 1.3.1 of
Robertson, Wright, and Dykstra, (1988, p.15);
and finally \begin{eqnarray*} \Sigma_2
& = & \sum_{g_*(s) > c} [g_*(s) - g_{cu}(s)][g_{cu}(s) -
		      \hat{g}(s)]w(s)\\
& = & \sum_{g_*(s) > c} [g_*(s) - c][c - \hat{g}(s)]w(s)
\geq 0 \;\;.  \end{eqnarray*} Since $\iuc$ is a convex
lattice we may apply the converse part of Theorem 1.3.1 of
Robertson, Wright, and Dykstra (1988),
and the result follows. \qued

Exactly analogous to Lemma 1 is

{\bf Lemma 2:}  The function \begin{displaymath}
g_{cl}(s) = \left \{ \begin{array}{cl} g_*(s) & \mbox{ if }
g_*(s) \geq c\\ c & \mbox{ if } g_*(s) < c \;\;.
\end{array} \right.  \end{displaymath} uniquely minimizes
(\ref{eq:trunciso}) for $\hat{g} \in \ilc$.

The following is an immediate consequence of Lemma 1 and 2:

{\bf Lemma 3:} Let $c_{k_1}, \ldots, c_{k_m}$ be a subset
of the level values of $g_*$, and let \linebreak $S' = S \setminus
\bigcup_{l=1}^m \{s | g_*(s) = c_{k_l} \} \neq \phi$.  The
isotonic regression of $g$ restricted to $S'$ is $g_*$
restricted to $S'$.

We can now prove the main result:

{\bf Proof of Theorem 1:} Since $S_1 \fol S_2$ it is easy to see that
there is a constant $c$ such that:  \begin{eqnarray*}
g_*(s) & < & c \mbox{ implies } s \in S_1 {\rm and}\\
g_*(s) & > & c \mbox{ implies } s \in S_2 \;\;.
\end{eqnarray*} The set $\{s | g(s) = c \}$ may contain
elements from both $S_1$ and $S_2$.

For this $c$ \begin{displaymath} g_*(s) = \left \{
\begin{array}{cl} g_{cu}(s) & \mbox{ if } s \in S_1 \\
g_{cl}(s) & \mbox{ if } s \in S_2  \end{array} \right.
\end{displaymath} otherwise we would contradict the
definition of $g_*$.  Applying lemmas 1 and 2, it follows
that \begin{displaymath} g_{cu}(s) = \left \{
\begin{array}{cl} g_{1*}(s) & \mbox{ if } g_{1*}(s) < c \\
c & \mbox{ if } g_{1*}(s) \geq c \end{array} \right.
\end{displaymath} for $s \in S_1$ and \begin{displaymath}
g_{lu}(s) = \left \{ \begin{array}{cl} c & \mbox{ if }
g_{2*}(s) \leq c \\ g_{2*}(s) & \mbox{ if } g_{2*}(s) > c
\end{array} \right.  \end{displaymath} for $s \in S_2$.
Therefore $g_*(s)$ is a function of $g_{1*}(s)$ on $S_1$,
and is a function of $g_{2*}(s)$ on $S_2$.  In other words,
$g_*(s)$ is constant on all of the level sets $L_i$ of
$g_{1*}$ and $g_{2*}$.  Let $g_*(s) = d_i$ on $L_i$.  Now
\begin{eqnarray*} \sum_S [g(s) - g_*(s)]^2w(s) & = &
\sum_{S_1} [g(s) - g_{1*}(s) + g_{1*}(s) - g_*(s)]^2w(s)\\
  & & + \sum_{S_2} [g(s) - g_{2*}(s) + g_{2*}(s) -
  g_*(s)]^2w(s)\\  & = & \sum_{S_1} [g(s) -g_{1*}(s)]^2w(s)
+ \sum_{S_2} [g(s) -g_{2*}(s)]^2w(s)\\
  & & + \sum_{S_1} [g_{1*}(s) -g_*(s)]^2w(s) + \sum_{S_2}
  [g_{2*}(s) -g_*(s)]^2w(s)\\ & & +  2 \sum_{S_1} [g(s) -
 g_{1*}(s)][g_{1*}(s) - g_*(s)]w(s)\\ & & + 2 \sum_{S_2}
 [g(s) - g_{2*}(s)][g_{2*}(s) - g_*(s)]w(s) \end{eqnarray*}
The last two terms are zero by Theorem 1.31 of
Robertson, Wright, and Dykstra, (1988) since
$g_{1*}(s) - g_*(s)$ is a function of $g_{1*}(s)$, and
$g_{2*}(s) - g_*(s)$ is a function of $g_{2*}(s)$.  The
first two terms do not involve $g_*(s)$.  Hence $g_*(s)$
minimizes \begin{equation} \sum_{S_1} [g_{1*}(s) -
g_*(s)]^2w(s) + \sum_{S_2} [g_{2*}(s) - g_*(s)]^2w(s)
\label{eq:minim} \end{equation} and hence is the isotonic
regression of \begin{displaymath} h(s) = \left \{
	\begin{array}{cl}
		g_{1*}(s) & \mbox{ if $s \in S_1$}\\
		g_{2*}(s) & \mbox{ if $s \in S_2$}\\
	\end{array} \right .  \end{displaymath} It follows
readily that the values of $g_*(s)$ on $L_i$, i.e. $d_i$,
are in increasing order.  Since $g_*(s)$ minimizes
(\ref{eq:minim}), equal to \begin{displaymath} \sum_{t=1}^{r}
[ c_t - d_t ]^2u(t) \end{displaymath} under the assumption
that $g_*$ is isotonic, it follows that $d_1, d_2, \ldots,
d_r$ minimize this expression under simple linear order on
$1, 2, \ldots, r$, and hence $d_t = f_*(t)$ for all $t$.
\qued

\begin{thebibliography}{99}

\bibitem{Gib} Gibbons, J. D., Olkin, I., and Sobel, M.
(1977).  {\em Selecting and Ordering Populations:  A New
Statistical Methodology.}  Wiley, New York.

\bibitem{Hot} Hotelling, H. (1941).  Experimental
determination of the maximum of a function.  {\em Ann. Math.
Statist.} {\bf 12}, 20-45.

\bibitem{Muk} Mukerjee, H.  (1988).  Monotone nonparametric
regression.  {\em Ann.  Statist.} {\bf 16}, 741-750.

\bibitem{Mul} M\"{u}ller, H.-G. (1989).  Adaptive
nonparametric peak estimation.  {\em Ann.  Statist.} {\bf 17},
1053-1069.

\bibitem{Par} Parzen, E. (1962).  On estimation of a
probability density function and mode.  {\em Ann. Math. Statist.}
{\bf 33}, 1065-1076.

\bibitem{Rob} Robertson, Tim.  (1978).  Testing for and
against an order restriction on multinomial parameters.  {\em J.
Amer. Statist. Assoc.} {\bf 73}, 197-202.

\bibitem{RWD} Robertson, T.,
Wright, F. T., and Dykstra, R.  L. (1988).  {\em Order Restricted
Statistical Inference.}  Wiley, New York.

\bibitem{Rom} Romano, J.  (1988).  On weak convergence and
optimality of kernel density estimates of the mode.  {\em Ann.
Statist.} {\bf 16} 629-647.

\bibitem{Sha} Shapiro, M.  (1985). Asymptotic distribution of
test statistics in the analysis of moment structures under
inequality constraints. {\em Biometrika,} {\bf 72} 133-144.

\bibitem{Sil} Silverman, B. W. (1985).  Some aspects of the
spline smoothing approach to nonparametric regression curve
fitting (with discussion).  {\em J. Roy. Statist. Soc. Ser. B}
{\bf 47}, 1-52.  \end{thebibliography}

\end{document}
