psMinimizeLMChi2()

[gusciora@…]

As a result of a conversation with Paul Price, we are no longer attempting a
general LMM minimization routine. Instead I have implemented a LMM chi-squared
minimization routine. Since function prototypes are not specified in the
current SDR, I implemented what seemed reasonable to me, however feel free to
modify it:

typedef
float (*psMinimizeLMChi2Func)(psVector *deriv,

const psVector *params,
const psVector *coords);

bool psMinimizeLMChi2(psMinimization *min,

psVector *params,
const psVector *paramMask,
psImage *covar,
const psArray *coords,
const psVector *value,
psMinimizeLMChi2Func func);

The significant difference is that psMinimizeLMChi2Func() takes as input a
single coord, rather than an array of coords. The reason for this is that the
algorithm in NR requires a derivative for each point, rather than a derivative
for all of the points, somehow. Let me know what you think.

[Paul Price]

typedef float (*psMinimizeLMChi2Func)(psVector *deriv, const psVector *params,
const psVector *x);

This function takes the current guess for the parameters for which we are trying
to get the best values (params), and a single vector (x) of conditions from the
array of x values fed into the minimiser.

bool psMinimizeLMChi2(psMinimization *min, psImage *covar, psVector *params,
const psVector *paramMask, const psArray *x, const psVector *y, const psVector
*yErr, psMinimizeLMChi2Func func);

This takes the minimization specs (min), returns the covariance matrix (covar),
takes the best guess of initial parameters (params), the parameter mask
(paramMask), and takes multiple vectors of conditions in an array (x), the
corresponding measured values (y) and errors (yErr) and the function to fit (func).

For example, for GRB afterglows, I have flux as a function of time and
frequency. So I stuff into the "psArray *x" all my time and frequency values
(so a whole heap of vectors of size 2), I have the measured values in "y" and
errors in "yErr". Then each of the time-frequency pairs are passed to my model
function with the current parameters, and the model function returns the flux
for that time-frequency pair, and the derivative with respect to each of the
parameters.

This seems reasonable. The only thing we could change would be to have the
function be defined:

typedef psVector* (*psMinimizeLMChi2Func)(psImage *deriv, const psVector
*params, const psArray *x);

So it would return the model value for each of the measurements at once, and
return for each the derivatives (so that it returns a matrix).

What do you think?

[Paul Price]

Updated the SDRS to reflect the above:

Consider a function of a collection of parameters, \code{params},
which is evaluated at a position, \code{x}, which returns a single
floating point value which is the value of the function given the
parameters and coordinate vectors, along with the derivatives of the
function with respect to each of the parameters, \code{deriv}:
\begin{verbatim}
typedef float (*psMinimizeLMChi2Func)(psVector *deriv, const psVector *params,
const psVector *x);
\end{verbatim}

Then \code{psMinimizeLMChi2} shall fit the specified function,
\code{func}, to a set of measurements, \code{x,y,yErr}, using the
Levenberg-Marquardt method:

\begin{verbatim}
bool psMinimizeLMChi2(psMinimization *min, psImage *covar, psVector *params,

const psVector *paramMask, const psArray *x, const

psVector *y,

const psVector *yErr, psMinimizeLMChi2Func func);

\end{verbatim}

The function shall return \code{false} in the event that there was an
error (bad input, too many iterations), and also call \code{psError}.

The minimization specification, \code{min}, shall be modified with the
\code{value}, \code{iter} and \code{lastDelta} members updated with
the values appropriate from the minimization.

On calling the minimizer, the \code{params} vector shall consist of a
best guess for the parameters that minimize the model function,
\code{func}. On successful completion, the input parameters,
\code{params}, shall be updated with the values that minimize the
input function. The function shall also update the covariance matrix,
\code{covar}, in the event that it is non-\code{NULL}. The function
shall generate in error in the event that \code{covar} is
non-\code{NULL} and is not a square matrix with size matching that of
\code{params}.

Parameters that have a corresponding \code{paramMask} entry of
\code{0} are to be held fixed by the minimizer. It shall be an error
for \code{paramMask} not to be of the same dimension as \code{params}.

The measurement ordinates, \code{x}, shall consist of multiple
vectors, each of which may be passed to the model \code{func}. If the
measurement coordinates, \code{y}, and errors, \code{yErr}, are not of
the same length as the ordinates array, \code{x}, then the function
shall generate a warning, and truncate the longest of the
array/vectors to match the length of the shortest. The vectors
contained within the \code{x} array, and the \code{y} and \code{yErr}
vectors must be of type \code{psF32}. The \code{yErr} vector may be
\code{NULL}, in which case the errors shall be assumed to be
identical.

\code{paramMask} must be of type \code{psU8}, while \code{params} must
be of type \code{psF32}. The \code{func} function must be valid only
for types \code{psF32}, \code{psF64}.

\subsubsubsection{Pre-defined Functions for LM}

We define some commonly used functions for use with the LM
minimization, used for the purpose of performing $\chi^{2$ fitting:}

\begin{verbatim}
psMinimizeLMChi2Func psMinimizeLMChi2Gauss1D;
psMinimizeLMChi2Func psMinimizeLMChi2Gauss2D;
\end{verbatim}

\code{psMinimizeChi2LMGauss1D} shall take as \code{params}, the
normalization, center, and standard deviation of a Gaussian to be fit,
and as \code{x}, a vector containing a single value. It shall return
the value of the Gaussian at the value, and the derivatives
(\code{deriv}) with respect to each of the parameters.

\code{psMinimizeChi2LMGauss2D} shall take, as \code{params}, the
normalization, center (two values), standard deviation (two values)
and position angle of a 2-dimensional Gaussian, and as \code{x}, a
vector containing a position, $(x,y)$. It shall return the value of
the 2-dimensional Gaussian at the specified point, along with the
derivatives with respect to each of the parameters.

[Paul Price]

Closing subsequent to release of SDRS-08, ADD-07.