\documentstyle[art10,npac_proto,palatino]{article}

\begin{document}
%\makeNPACtitle\newpage
\begin{project}
\begin{titlepage}
	\projtitle{NPAC Sailfish}
	\doctitle{User Guide for the Stochastic Global Optimisation
Package `Laplace'}
		\author{Miloje S. Makivic}
		\date{August 21, 1995}
		\version{0.2}
		\summary{Delivered August 23, 1995}
	\id{laplace}
	\docauthor{M~S~Makivic}
%	\doccontributor{KAH}
	\security{Commercial in Confidence}
\maketitle
	\approval{NPAC}
\vfill
	\begin{history}
		\release{
			\author{Miloje S. Makivic}
			\date{August 3, 1995}
			\version{0.1}
			\summary{Initial draft}
			}
	\end{history}
\end{titlepage}

\pagenumbering{roman}\setcounter{page}{1}
\begin{abstract}
`Laplace' is a stochastic global optimization package based on the classical Laplace formula.
The non-degenerate maximum or minimum value of a function is computed as a 
limit $\lambda\to\infty$ of the
first moment of the exponent of the objective function multiplied by the factor $\lambda$.
In addition to the location of extremum, the package computes error bounds 
for maximum likelihood or least squares estimation. 
\end{abstract}

{\parskip=0pt   % hack to patch format of \tableofcontents
\tableofcontents
%\newpage
\medskip\hspace{2.0cm}\hrulefill\hspace{2.0cm}\bigskip
}

\pagenumbering{arabic}\setcounter{page}{1}

%=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
\section{Introduction}
One of the principal methods of time series analysis is the Maximum Likelihood Estimation
of process parameters [1]. 
Model of stochastic rate movement involves either a probability distribution
for noise or an underlying Markov process. To estimate the parameters of this probability
distribution (we will refer to this task as calibration), 
one first constructs the log-likelihood function using an observed sample
of process history and the postulated process probability distribution. 
The set of parameters of the probability distribution which maximizes the log-likelihood
function is the maximum likelihood parameter estimate.
Therefore, the central task in maximum likelihood estimation is the maximization of log-likelihood
function.
The topology of the log-likelihood function can be quite complex, so robust methods for large scale
global optimization must be used. Global optimization methods are required, since there may be
multiple local extrema, depending on the function itself and process sample size.
NPAC provides two global stochastic optimization packages: one is based on Laplace formula and is
similar in spirit to simulated annealing, while the other one uses a stochastic differential
equation solver. This document describes the calibration package based on  Laplace formula. 
\section{Description of Algorithm}
	The algorithm [2] finds the global
maximum of a function $f(\vec{x})$ of
$\vec{x}=(x_1,x_2,...,x_i ,...,x_n)$ degrees of freedom. 
This function will be referred to as the objective function in the remainder of this guide.
Location of the maximum is given by the Laplace formula:
\begin{equation}
\vec{x}^* = \lim_{\lambda\to\infty} 
               {
               {\int_D \vec{x}\exp(\lambda f(\vec{x}))d\vec{x}} 
                       \over 
               {\int_D \exp(\lambda f(\vec{x}))d\vec{x}}
               } = \lim_{\lambda\to\infty} <\vec{x}>
\label{laplace}
\end{equation}
The estimation of the first moment $<\vec{x}>$ for each multiplier $\lambda$
is done using the Metropolis Monte Carlo
algorithm. Simulation begins with multiplier $\lambda = 1$ and proceeds
for larger values of multipliers as specified by the user. 
The estimate $<\vec{x}>$ from previous
multiplier is used as the initial point for the simulation with next
multiplier. 

By increasing $\lambda$, vectors $\vec{x}$ generated
in the Metropolis algorithm are more tightly bound to the location of maximum.
Their spread around the average value is a measure of the uncertainty
of the location of maximum. This uncertainty is reduced to $0$ when $\lambda\to\infty$.
Since this limit cannot be achieved in practice, 
user has to specify the desired accuracy of the maximum location.
Two convergence criteria are available:
\begin{itemize}
\item strong convergence
\item weak convergence
\end{itemize}

If the user requests strong convergence criterion,
simulation terminates when the Monte Carlo estimate of
variance for each degree of freedom, $\sigma_i = (<x_i^2 > - <x_i>^2)^{1/2}$ becomes smaller than user
specified accuracy. The accuracy is specified as a small parameter $\epsilon_i$ for each 
degree of freedom $x_i$ in the input
file. The criterion has different forms for degrees of freedom above and
below a user specified threshold $x_i^T$:
\begin{itemize}
\item Strong relative convergence: $\sigma_i \le \epsilon_i \cdot <x_i>$ if
$<x_i> \le x_i^T$
\item Strong absolute convergence: $\sigma_i \le \epsilon_i$ if $<x_i> > x_i^T$
\end{itemize}
These different criteria ensure that time is not wasted trying to converge for near zero parameters,
which are likely to be irrelevant anyway.
The user should ensure that the largest multiplier value $\lambda_{max}$
used in the optimization process is sufficient to obtain the desired accuracy levels.
Otherwise, simulation stops at $\lambda_{max}$ and its accuracy is determined by $\lambda_{max}$.

The weak convergence option for accuracy control is useful for functions with
either high degree of symmetry near the extremum or flat extrema.
If the difference between $\vec{x}^* = <\vec{x}>$ estimates for two 
successive multipliers $\lambda_1 < \lambda_2$ is smaller than the user specified threshold,
the algorithm is considered to have weakly converged. More precisely, user specifies desired
accuracy as an input parameter $\epsilon$.
If $| <x_i (\lambda_2 ) > - <x_i (\lambda_1 )> | \le \epsilon \sigma_i ((\lambda_2 )$, for each
degree of freedom $x_i$, 
weak convergence is achieved.
 It is user's responsibility to
ensure that $\lambda_2$ is sufficiently bigger than $\lambda_1$, otherwise
this convergence criterion is meaningless. In general, this criterion should be
used only if user has some knowledge about function behavior near extremum. 

An important parameter for the Metropolis algorithm is the typical size of a
perturbation (trial update) around the current configuration. For the initial
simulation at $\lambda=1$, this update range is specified by the user. For
all other values of $\lambda$, it is determined using the variance $\sigma_i$ as computed
for the previous value of $\lambda$. This ensures that ranges of trial updates for each
degree of freedom reflect the
scale of variation of the objective function near its maximum. 

In the context of time series analysis, the Laplace package will be used
for maximum likelihood estimation, although it can be used for any kind of large
scale optimization problem. 
User has an option to mask
objective function function degrees of freedom which are not variational. This mask is processed
by the Metropolis update routine which ensures that trial moves
are not attempted along the masked coordinates. 

The package provides a sample interface to
presently used NPAC/Sailfish data structures for time series modeling. Also,
the interface to the objective function contains a pointer to the PATH
structure with time series data. For purposes of time series estimation,
Laplace package computes parameter error estimates using standard procedures.
Log-likelihood function is regarded as a probability function with
Gaussian behavior near the maximum. Parameter error bounds are then obtained
as diagonal elements of the covariance matrix, $(<x_i^2 > - <x_i>^2)^{1/2}$.
Monte Carlo can easily be extended to compute full error covariance matrix,
$(<x_i x_j > - <x_i><x_j >)^{1/2}$, but that would incur quite high
CPU and memory costs.

%=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
\section{Description of Code}
	At the time of writing, the Laplace package has been compiled
and tested for the following architectures: Sparc Sun~4, IBM RS/6000
AIX~3.2, and Digital Alpha OSF/1~V3.0.  Compilation of the package is
via the standard Unix {\tt make} command. 
The code is
written in ANSI-C, and compiled with the GNU gcc compiler.

Laplace package consists of following modules:
\begin{itemize}
	\item {\tt modest.c}, sample routine to issue requests and obtain
reports for time series maximum likelihood estimation. This module
relies on services provided by the optimizer routines in {\tt modlaplace.c}.

	\item {\tt modlaplace.c}, routines for finding global extremum of
a function using Laplace formula and Metropolis Monte Carlo simulation.
To access services of the Laplace optimizer one first issues a request
to create an optimizer. The client application receives a handle to the
optimizer which is later used to retreive information from the optimizer.
Multiple optimizers, differentiated by their unique hadles, can coexist.

	\item {\tt modmc.c}, routines for Metropolis Monte Carlo simulation. 
        \item {\tt moments.c}, Routines for accumulating first and second
moment during MC simulation.
	\item {\tt modrng.c}, routines for random number generation. 
	\item {\tt modpath.c}, data structures and routines to manipulate
time series data. This is the same module as the one used by the Monte Carlo Simulator
for proxy history generation.

	\item {\tt memory.c}, utility routines for memory management.
	\item {\tt modvec.c}, utility routines for manipulation of vectors and matrices.
	\item {\tt modproc.c}, data structures and routines for manipulating
parameters of the stochastic process. This module is a generaliztion of the {\tt modproc.c} 
module used by the Monte Carlo Simulator for proxy history generation. 

	\item {\tt modfunc.c}, user supplied library of objective functions to be minimized.
Log-likelihood functions for time series parameter estimation should be defined
in this module.

	\item {\tt modutil.c}, miscellaneous utility functions.

\end{itemize}
In the following sections, functionality of important modules will be
described in detail.

%=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
\section{Module MODEST}

The main routine in this module is called {\tt EstimateProcessParameters}.
It illustrates how to use services of the optimizer in order to obtain
estimates of the process parameters and estimation errors and generate
estimation reports.
The calling sequence is:
\begin{verbatim}
int   EstimateProcessParams
                          (
                          PATH *Noise,                  /* input  */
                          PROCESS_EST *ProcessParams,   /* in/out */
                          char *Optimizer_descriptor    /* input  */
                          )

\end{verbatim}
It takes as input a sample of process history ({\bf Noise}), a data structure which
describes process parameters ({\bf ProcessParams}) and the name of the
file which describes parameters to configure the optimizer ({\bf Optimizer\_descriptor}).
{\bf ProcessParams} contains initial guesses for parameter estimates and their
uncertainties, the address of
the maximum likelihood function, and it also
identifies which parameter values should be fixed and which ones are variational.
Upon completion, the appropriate  fields in {\bf ProcessParams} will contain 
maximum likelihood estimate of parameters and
their estimated errors.

These are the basic steps in {\tt EstimateProcessParams} (see also comments in the source code): 
\begin{itemize}
	\item Assemble input variables for the optimizer from {\bf ProcessParams} 
        \item Create an optimizer
        \item Invoke the optimizer to maximize likelihood function
        \item Obtain parameter error estimates from the optimizer
        \item Generate report
        \item Destroy the optimizer
\end{itemize}

	The user has to provide definition of log-likelihood function (which is the objective
function in this context) with the
following interface:

\begin{verbatim}
   (REAL) LogLikelihood
                      (
                      int NumOfArguments, 
                      REAL *ArgumentArray, 
                      PATH *ProcessSample
                      )
\end{verbatim}

where the C-preprocessor constant `{\tt REAL}' specifies
whether type {\tt double} or {\tt float} is used; `{\tt NumOfArguments}' is the
number of degrees of freedom of the log-likelihood function; and `{\tt ArgumentArray}'
is a pointer to a {\tt NumOfArguments}-dimensional vector of the degrees of
freedom. The optimizer can be used as is for any kind of optimization problem by
setting the PATH pointer variable to NULL.

  
\section{Module MODLAPLACE}

The services of MODLAPLACE module are generally accessed through the following sequence of functions:
\begin{itemize}
\item Creation of optimizer,
\begin{verbatim}

       int CreateOptimizer
                         (
                         int NumArgs,
                         char *Optimizer_descriptor
                         )
\end{verbatim}
{\tt NumArgs} is the number of degress of freedom for optimization, and {\tt Optimizer\_descriptor}
is the name of the file which is used to configure the optimizer. The description of input
parameters and input file formats is given in the Laplace Optimizer Initialization Section. {\tt CreateOptimizer} returns
an integer handle to an OPT\_LAPLACE data structure (described in the Data Strcuctures section). 
Each call to {\tt CreateOptimizer} allocates a new optimizer which is accessed by the client
application through its
unique integer handle.
\item Optimization request,
\begin{verbatim}
        int  MaximizeObjectiveFunction
                                     (
                                     int OptimizerID,   
                                     PATH *Noise,
                                     int NumArgs,
                                     REAL *Args,
                                     REAL *Sigmas,
                                     int *VarArgs,
                                     ObjFunc *ObjectiveFunction
                                     ) 
\end{verbatim}
{\tt OptimizerID} is the optimizer handle returned by a previous call to {\tt CreateOptimizer};
{\tt Noise} is a sample of process history; {\tt NumArgs} is the number of degress of freedom for optimization
(same as the one specified to {\tt CreateOptimizer}). {\tt Args} is the array which stores the
initial guess for location of the maximum when {\tt MaximizeObjectiveFunction} is called and
the final result for location of the maximum upon return from {\tt MaximizeObjectiveFunction}.
In the context of maximum likelihood estimation this would be the initial guess for process
parameters.
{\tt Sigmas} contains the initial guess for errors of process parameter estimates. It is not
necessary to have a good guess for either {\tt Args} or {\tt Sigmas}. After the initial
simulation with multiplier $\lambda=1$, both of these are automatically adjusted by the
optimizer. The initial guess for {\tt Sigmas} is used to set the initial scale for Metropolis
updates. Later on it is adjusted as described in the algorithm Section . {\tt VarArgs} is the
mask which specifies whether a degree of freedom is variational or not: if {\tt VarArgs[i] = 1}
it is variational, otherwise it is fixed to its initial value in the course of the optimization. 
This is the only constraint for optimization which can be specified as an argument to the
{\tt MaximizeObjectiveFunction} routine. Any other constraint must be specified through the
objective function itself, which is the last argument, {\tt ObjectiveFunction}. In the
context of time series estimation, this would be the log-likelihood function.
\item Computation of  estimation errors,
\begin{verbatim}
      int   MaxLikelihoodErrorEstimate
                                     (
                                     int OptimizerID,
                                     int NumArgs,
                                     REAL *Sigmas
                                     )         
\end{verbatim}
{\tt OptimizerID} is the optimizer handle returned by a previous call to {\tt CreateOptimizer};
{\tt NumArgs} is the number of degress of freedom for optimization 
(same as the one specified to {\tt CreateOptimizer} and {\tt MaximizeObjectiveFunction}).
The error estimate is returned in the {\tt Sigmas} array.
\item Report generation,
\begin{verbatim}
      int   EstimationReport
                           (
                           int OptimizerID
                           )
\end{verbatim}
{\tt OptimizerID} is the optimizer handle returned by a previous call to {\tt CreateOptimizer}.
The optimizer data structure stores all the information relevant for the analysis of the
optimization procedure, including the type of report requested in the input file. The
presently available choices are:
\begin{enumerate}
\item No report at all
\item Short report, which gives only location of maximum and its uncertainty
\item Long report, which gives additional information about convergence and consistency of
     optimization
\end{enumerate}
\item Destruction of optimizer,
\begin{verbatim}
      void    CloseOptimizer
                           (
                           int OptimizerID
                           )
\end{verbatim}
\end{itemize}


\section{Module MODMC}
This is the Metropolis Monte Carlo simulation module which computes $<\vec{x}>$ and 
$<\vec{x}^2 > - <\vec{x}>^2$ for each multiplier $\lambda$. 
The interface to the module is exactly the same as for the Metropolis Monte Carlo Simulator 
for proxy rate generation (Proxy Rate Simulator), so the reader is referred to the user guide provided earlier
for the Proxy Rate Simulator.
Since a Metropolis simulation must be performed for each value of the multiplier $\lambda$,
user has to specify the Metropolis driver descriptor file to be used for each value of $\lambda$.
The file names will be stored in the optimizer OPT\_LAPLACE data structure and accessed as needed
during the optimization procedure. It is possible to use the same input file for all simulations.
However, it is likely that a user will have a very poor initial guess for the maximum location.
It may also happen that the initial, user-specified scale for updates results in a poor efficiency of
the Metropolis algorithm.
Therefore, it is recommended that the simulation for $\lambda=1$ is run with a bigger number
of Monte Carlo steps than other multipliers, which are going to inherit valuable information
about the objective function from the $\lambda=1$ simulation.

The simulation uses the same basic
algorithm as the Metropolis Monte Carlo Simulator for proxy rate generation, but the
implementation is rather different to reflect the different needs of the optimizer.
The objective function is evaluated for initial and final configuration of the
Metropolis update, unlike the Proxy Rate Simulator, which evaluates the transition probability
directly given the initial configuration and its perturbation. 
It is also possible to mask the coordinates which should not be updated. Finally,
the Laplace driver for Metropolis update algorithm is now capable of 
automatic rescaling of the range of updates for each
$\lambda$, since the user is not expected to be able to provide this
parameter accurately. Users can pick between 2 alternatives for automatic
rescaling using the pre-processor {\tt DEFINE} variable {\tt RESCALE}:
\begin{itemize}
\item If {\tt RESCALE} is {\tt DEFINED}, the range of updates for a
current multiplier value is computed as the variance from previous multiplier
rescaled by the square root of the ratio of $\lambda$ multipliers. This will
typically result in higher acceptance ratios and more correlated configurations.
It is preferred for sharp increase between successive $\lambda$ multipliers.
\item If {\tt RESCALE} is not {\tt DEFINED}, the range of updates for a
current multiplier value is computed as the variance from previous multiplier.
This will typically result in a smaller acceptance ratio but less
correlated configurations. This is the preferred method for slower
rate of increase of $\lambda$ multipliers.  
\end{itemize}

\section{Data Structures}

The central data structure of the calibration module is {\tt OPT\_LAPLACE}:
\begin{verbatim}
typedef struct opt_laplace{
        int Final;
        int TypeOfReport;
        int Dimension;
        int NumMultipliers;
        REAL UpdateRescaleFac;
        int RequestStrongConvergence;
        int StrongConvergence;
        int WeakConvergence;
        REAL WeakConvergenceEpsilon;
        REAL ArgThreshold;
        REAL AbsEpsilon;
        REAL RelEpsilon;
        REAL MaxFuncVal;
        REAL MinFuncVal;
        REAL Maximum;
        int CurrentMX;
        int LastUsedMX;
        REAL *Multipliers;
        REAL MaxMultiplier;
        char **MC_Descriptors;
        char *ReportFile;
        REAL **Avg;
        REAL **Sigma;
} OPT_LAPLACE;

\end{verbatim}
The fields of this data structure are defined as follows:
\begin{itemize}
\item {\bf Final}, a flag which specifies whether the optimizer is used to generate final
result or as a preprocessor to generate initial starting point for another optimizer.
In the latter case, convergence test is not performed.
\item {\bf TypeOfReport}, this flag identifies the type of report requested from the
optimizer. if set to $0$, it will not produce report at all, if set to $1$ it will produce
a short report, and if set to $2$ it will produce a detailed report.
\item {\bf Dimension}, number of degrees of freedom in the objective function,
\item {\bf NumMultipliers}, number of multipliers $\lambda$, in addition to the unit
multiplier, to be used in the optimization,
\item {\bf UpdateRescaleFac}. The user is required to provide an estimate of uncertainty
for the initial guess of optimal parameters. This factor is used to rescale this guess
if deemed necessary.
\item {\bf RequestStrongConvergence}, an integer flag which specifies whether strong
convergence must be achieved for successful completion. Weak convergence criterion
is the default, so if this flag is set to $0$, optimizer will exit as soon as weak
convergence is achieved. If this flag is set to $1$, both weak and strong convergence
criteria will be tested, but successful completion requires that strong criterium is met.
\item {\bf StrongConvergence}, this flag indicates whether strong convergence is achieved,
\item {\bf WeakConvergence}, this flag indicates whether weak convergence is achieved,
\item {\bf WeakConvergenceEpsilon}, this number defines the accuracy required by the
weak convergence criterion. 
\item {\bf ArgThreshold}, for degrees of freedom which are below this value, absolute
strong convergence is requested, while relative strong convergence applies for those above
threshold.
\item {\bf AbsEpsilon}, absolute accuracy required for absolute strong convergence criterion.
\item {\bf RelEpsilon}, relative accuracy required for relative strong convergence criterion.
\item {\bf MaxFuncVal}, maximum value of the objective function encountered in the course of
Monte Carlo simulation.
\item {\bf MinFuncVal}, minimum value of the objective function encountered in the course of
Monte Carlo simulation.
\item {\bf Maximum}, value of objective function at its estimated maximum.
\item {\bf CurrentMX}, index of the current value of the multiplier $\lambda$
\item {\bf LastUsedMX}, index of the last multiplier used in the course of optimization.
\item {\bf Multipliers}, the array of $\lambda$ multiplier values to be used in optimization.
\item {\bf MaxMultiplier}, the biggest multiplier used in optimization. If no convergence
was accomplished it is equal to the biggest multiplier requested by the user.
\item {\bf MC\_Descriptors}, list of files with input parameters to configure Metropolis drivers.
One file name must be specified for each multiplier, although they do not have to be all different.
\item {\bf ReportFile}, name of the report file.
\item {\bf Avg}, storage for array $<\vec{x} >$ for each multiplier.
\item {\bf Sigma}, storage for array $\vec{\sigma} = <\vec{x}^2 > -<\vec{x}^2 >^{1/2}$ for each multiplier.

\end{itemize}

The interface between the calibration module and simulation modules is
achieved via the {\tt PROCESS\_EST} data structure which is derived from {\tt PROCESS} data
structure used in the Monte Carlo Simulator. {\tt PROCESS\_EST} has a couple of extra fields
to record errors for parameter estimates and mask for specification of variational parameters.
Also, instead of a pointer to the probability density function which defines the process, it
stores a ponter to the maximum likelihood function. Log-likelihood function
$\Omega (\vec{\theta})$ and process
probability density $f(\vec{x}_t , \vec{\theta} )$ 
are related via equation:
$$
\Omega(\vec{\theta}) = \sum_{t=1}^{t=T} \ln f(\vec{x}_t , \vec{\theta} )
$$
where $\vec{x}_t$ is the observed history at time step $t$ and $\vec{\theta}$
is the sought parameter array.
\begin{verbatim}
typedef struct process_est {
       int ProcessID;
       int Dimension;
       ProbFunc *Probability;
       int *CouplingsMap;
       REAL **RateCouplings;
       REAL **CouplingErrors;
       int  **VarCouplings;
} PROCESS_EST;
\end{verbatim}
The fields of this data structure are defined as follows:
\begin{itemize}
\item {\bf ProcessID}, process index which is used to identify the process likelihood function. 
\item {\bf Dimension}, number of instruments (rates) which are observed in the process.
\item {\bf Likelihood}, log-likelihood probability function.
\item {\bf CouplingsMap}, a list of number of couplings for each instrument (rate). For a detailed
explanation, consult the user guide to data structures of the Proxy Rate Simulator.
\item {\bf RateCouplings}, storage for instrument couplingss. These arrays define the
parameters of both process probability density and log-likelihood function. For a detailed
explanation, consult the user guide to data structures of the Proxy Rate Simulator.
\item {\bf CouplingErrors}, storage for estimated errors of instrument couplings. There is one-to-one
correspondenec between elements of these arrays and {\tt RateCouplings} arrays.
\item {\bf VarCouplings}, storage for flags indicatng whether a parameter is variational (flag = $1$)
or fixed (flag = $0$). There is one-to-one
correspondenec between elements of these arrays and {\tt RateCouplings} arrays.
\end{itemize}

The remaining major data structure is {\tt PATH}, which is the same as in
the Proxy Rate Simulator. The {\tt MODPATH} modules for calibration and simulation are also identical.



\section{Laplace Optimizer Initialization}
Input parameters for creation and initialization of the optimizer are given in the input file
with following entries:
\begin{itemize}
\item The {\bf Final} flag
\item Type of report {\bf TypeOfReport}
\item Report file name, {\bf ReportFile}
\item {\bf ArgThreshold}
\item {\bf AbsEpsilon}
\item {\bf RelEpsilon}
\item {\bf UpdateRescaleFactor}
\item flag for strong convergence, {\bf StrongConvergenceRequested}
\item {\bf WeakConvergenceEpsilon}
\item number of non-unit exponential multipliers, {\bf NumMultipliers} 
\item a list of multiplier values, followed by a filename with input parameters
for the Metropolis simulation with that multiplier value.
\end{itemize}

Format of the Metropolis simulation input files is described in the simulator documentation.

\section{Sample Program}
The distribution directory for the Laplace optimizer package contains an example program
which maximizes an inverted quadratic function. The example code is compiled by typing
{\tt make sail}. To run the program, one should type {\tt sail test.in test.out opt.in}.
Three input files have to be prepared: 
\begin{itemize}
\item {\tt opt.in}: optimizer initializer file with entries as described in previous Section
\item  {\tt mc.in}: Metropolis initializer file with entries as described in the Simulator
documentation; for the example code, simulations for all multipliers use this same input file.
This file name is specified in the {\tt opt.in} file.
\item {\tt test.in}: this input file serves as a proxy for real process parameters. In a real
application, the {\tt PROCESS\_EST} structure would be passed to the {\tt EstimateProcessParams}
routine. {\tt EstimateProcessParams} would in turn invoke the Laplace optimizer.
In the example code, the {\tt test.in} file is processed by the {\tt TestOptimizer}
routine in {\tt modest.c}, which invokes the Laplace optimizer to find the maximum of a 
quadratic function centered at $\vec{x} = 0$.
The format of {\tt test.in} is self-explanatory. 
\end{itemize}

The output of example code is written to {\tt test.out}. Also, if a report is requested by
specifying the appropriate flag in {\tt opt.in}, the report file {\tt report.out} is generated.
File name {\tt report.out} is one of the parameters specified in {\tt opt.in}.

Contents of {\tt report.out} depend on the type of report requested. For the
detailed report, it will specify what type of convergence is achieved. Then it will give
the largest value of the multiplier $\lambda$ used in the simulation. Following that,
it will list maximum and minumum value as recorded in the Monte Carlo simulation. These
are not statistical estimates. The Metropolis routine computes objective function
value for each configuration reached during the simulation and it will record the
max and min values obtained. This is a good consistency check for the optimization, since
the estimated maximum value must be equal or smaller than any of the values obtained in the
simulation. Otherwise, the simulation did not relax to the global maximum and should be
repeated with a different random number seed. 
The estimated maximum value is printed immediately after the recorded max and
min values. Following that is a table which lists components of the maximum vector as
estimated for each $\lambda$ value used in the simulation. The last two columns have
the following headings: Location Error and Estimation Error. Location error is given
by the variance of $\vec{x}$ at a given $\lambda$. Estimation Error is obtained by
multiplying Location Error with $\lambda^{1/2}$ and it is equal to the diagonal
element of the inverse curvature matrix of the objective function at the estimated
maximum. Estimation error is also recorded in the {\tt PROCESS\_EST} structure, along
with the location of maximum.
 

\section{References}
\begin{itemize}
\item[[1]] H. Lutkepohl, 1991, {\it Introduction to Multiple Time Series Analysis} (Springer-Verlag).
\item[[2]] Pincus, M., 1970, A Monte Carlo Method for the Approximate Solution of Certain Types of
Constrained Optimization Problems, {\it Operations Research}, {\bf 18}, 1225-1228.

\end{itemize}



\end{project}
\end{document}


