Fits a multivarite nonlinear regression model.
Synopsis
#include <imsls.h>
float *imsls_f_nonlinear_regression (float fcn(), int n_parameters, int n_observations, int n_independent, float x[], float y[], ..., 0)
The type double function is imsls_d_nonlinear_regression.
Required Arguments
float fcn (int n_independent, float xi[], int n_parameters,
float theta[])
User-supplied
function to evaluate the function that defines the nonlinear regression problem
where xi is an
array of length n_independent at which
point the function is evaluated and theta is an array of
length n_parameters
containing the current values of the regression coefficients. Function fcn returns a
predicted value at the point xi. In the following,
f(xi;θ), or just
fi, denotes the value of this function at the point
xi, for a given value of θ. (Both xi
and θ are arrays.)
int
n_parameters (Input)
Number of parameters to be
estimated.
int
n_observations (Input)
Number of observations.
int n_independent
(Input)
Number of independent variables.
float x[]
(Input)
Array of size n_observations by
n_independent
containing the matrix of independent (explanatory) variables.
float y[]
(Input)
Array of length n_observations
containing the dependent (response) variable.
Return Value
A pointer to an array of length n_parameters
containing a solution, 
for the nonlinear regression coefficients. To release this
space, use imsls_free.
If no solution can be computed, then NULL
is returned.
Synopsis with Optional Arguments
#include <imsls.h>
float
*imsls_f_nonlinear_regression (float
fcn(),
int n_parameters, int n_observations, int
n_independent, float x[], float
y[],
IMSLS_THETA_GUESS, float
theta_guess[],
IMSLS_JACOBIAN, void
jacobian(),
IMSLS_THETA_SCALE, float
theta_scale[],
IMSLS_GRADIENT_EPS, float
gradient_eps,
IMSLS_STEP_EPS, float step_eps,
IMSLS_SSE_REL_EPS, float sse_rel_eps,
IMSLS_SSE_ABS_EPS, float sse_abs_eps,
IMSLS_MAX_STEP, float max_step,
IMSLS_INITIAL_TRUST_REGION, float
trust_region,
IMSLS_GOOD_DIGIT, int ndigit,
IMSLS_MAX_ITERATIONS, int max_itn,
IMSLS_MAX_SSE_EVALUATIONS, int max_sse_eval,
IMSLS_MAX_JACOBIAN_EVALUATIONS, int max_jacobian,
IMSLS_TOLERANCE, float tolerance,
IMSLS_PREDICTED, float **predicted,
IMSLS_PREDICTED_USER, float predicted[],
IMSLS_RESIDUAL, float **residual,
IMSLS_RESIDUAL_USER, float residual[],
IMSLS_R, float **r,
IMSLS_R_USER, float r[],
IMSLS_R_COL_DIM, int r_col_dim,
IMSLS_R_RANK, int *rank,
IMSLS_X_COL_DIM, int x_col_dim,
IMSLS_DF, int *df,
IMSLS_SSE, float *sse,
IMSLS_VARIANCE_COVARIANCE_MATRIX, float **var_covar,
IMSLS_VARIANCE_COVARIANCE_MATRIX_USER, float var_covar[],
IMSLS_RETURN_USER, float theta_hat[],
IMSLS_FCN_W_DATA, void fcn(),void
*data,
IMSLS_JACOBIAN_W_DATA, void jacobian(),
void
*data,
0)
Optional Arguments
IMSLS_THETA_GUESS, float
theta_guess[] (Input)
Array with n_parameters
components containing an initial guess.
Default: theta_guess[] = 0.
IMSLS_JACOBIAN, void jacobian
(int
n_independent,
float
xi[], int n_parameters,
float theta[], float fjac[])
(Input/Output)
User-supplied function to compute the i-th row of the
Jacobian, where the n_independent data
values corresponding to the i-th row are input in xi. Argument theta is an array of
length n_parameters
containing the regression coefficients for which the Jacobian is evaluated,
fjac is the
computed n_parameters row of
the Jacobian for observation i at theta. Note that each
derivative ∂f(xi)/∂θj should be
returned in fjac
[j − 1] for
j = 1, 2, ..., n_parameters.
IMSLS_THETA_SCALE, float
theta_scale[] (Input)
Array with n_parameters
components containing the scaling array for θ. Array theta_scale is used
mainly in scaling the gradient and the distance between two points. See keywords
IMSLS_GRADIENT_EPS and
IMSLS_STEP_EPS
for more details.
Default: theta_scale[] = 1.
IMSLS_GRADIENT_EPS, float
gradient_eps (Input)
Scaled gradient tolerance. The
j-th component of the scaled gradient at θ is calculated as
where g = ∇F(θ), t = theta_scale, and
The value F(θ) is
the sum of the squared residuals, SSE, at the point θ.
Default:
(
in double, where
ɛ is the machine
precision)
IMSLS_STEP_EPS, float step_eps
(Input)
Scaled step tolerance. The j-th component of the scaled step
from points θ and θʹ is computed as
where t = theta_scale
Default:
step_eps = ɛ2/3,where ɛ is the machine precision
IMSLS_SSE_REL_EPS, float
sse_rel_eps (Input)
Relative SSE function
tolerance.
Default: sse_rel_eps = max(10-10, ɛ2/3), max(10-20, ɛ2/3) in double, where ɛ is the machine precision
IMSLS_SSE_ABS_EPS, float sse_abs_eps
(Input)
Absolute SSE function tolerance.
Default: sse_abs_eps = max(10-20,ɛ2), max(10-40, ɛ2) in double, where ɛ is the machine precision
IMSLS_MAX_STEP, float max_step
(Input)
Maximum allowable step size.
Default: max_step = 1000 max
(ɛ1, ɛ2), where ɛ1 = (tTθ0)1/2, ɛ2 = ||t||2,
t = theta_scale, and θ0 = theta_guess
IMSLS_INITIAL_TRUST_REGION, float
trust_region (Input)
Size of initial trust region radius.
The default is based on the initial scaled Cauchy step.
IMSLS_GOOD_DIGIT, int ndigit
(Input)
Number of good digits in the function.
Default: machine
dependent
IMSLS_MAX_ITERATIONS, int max_itn
(Input)
Maximum number of iterations.
Default: max_itn = 100
IMSLS_MAX_SSE_EVALUATIONS, int
max_sse_eval (Input)
Maximum number of SSE function
evaluations.
Default: max_sse_eval = 400
IMSLS_MAX_JACOBIAN_EVALUATIONS, int
max_jacobian (Input)
Maximum number of Jacobian
evaluations.
Default: max_jacobian = 400
IMSLS_TOLERANCE, float tolerance
(Input)
False convergence tolerance.
Default: tolerance = 100* eps, where
eps = imsls_f_machine(4) if
single precision and eps = imsls_d_machine(4) if
double precision
IMSLS_PREDICTED, float
**predicted (Output)
Address of a pointer to a real
internally allocated array of length n_observations
containing the predicted values at the approximate solution.
IMSLS_PREDICTED_USER, float
predicted[] (Output)
Storage for array predicted is provided
by the user. See IMSLS_PREDICTED.
IMSLS_RESIDUAL, float
**residual (Output)
Address of a pointer to a real
internally allocated array of length n_observations
containing the residuals at the approximate solution.
IMSLS_RESIDUAL_USER, float
residual[] (Output)
Storage for array residual is provided
by the user. See IMSLS_RESIDUAL.
IMSLS_R, float **r
(Output)
Address of a pointer to an internally allocated array of size
n_parameters × n_parameters
containing the R matrix from a QR decomposition of the
Jacobian.
IMSLS_R_USER, float r[]
(Output)
Storage for array r is provided by the
user. See IMSLS_R.
IMSLS_R_COL_DIM, int r_col_dim
(Input)
Column dimension of array r.
Default: r_col_dim = n_parameters
IMSLS_R_RANK, int *rank
(Output)
Rank of r. Argument rank less than n_parameters may
indicate the model is overparameterized.
IMSLS_X_COL_DIM, int x_col_dim
(Input)
Column dimension of x.
Default: x_col_dim = n_independent
IMSLS_DF, int *df
(Output)
Degrees of freedom.
IMSLS_SSE, float *sse
(Output)
Residual sum of squares.
IMSLS_RETURN_USER, float
theta_hat[] (Output)
User-allocated array of length n_parameters
containing the estimated regression coefficients.
IMSLS_VARIANCE_COVARIANCE_MATRIX,
float
**var_covar (Output)
Address of a pointer to an internally
allocated array of size n_parameters × n_parameters
containing the variance/covariance matrix of the coefficients.
IMSLS_VARIANCE_COVARIANCE_MATRIX_USER,
float
var_covar[] (Output)
Storage for array var_covar is provided
by the user. See IMSLS_VARIANCE_COVARIANCE_MATRIX.
IMSLS_FCN_W_DATA, float fcn (int n_independent, float xi[], int n_parameters,
float theta[]),
void
*data, (Input)
User-supplied function to evaluate the function that
defines the nonlinear regression problem, which also accepts a pointer to data
that is supplied by the user. data is a pointer to
the data to be passed to the user-supplied function. See the “Introduction, Passing Data to
User-Supplied Functions” at the beginning of this manual for more details.
IMSLS_JACOBIAN_W_DATA, void jacobian
(int
n_independent,
float
xi[], int n_parameters,
float theta[], float fjac[]),
void
*data, (Input)
User-supplied function to compute the i-th row
of the Jacobian, which also accepts a pointer to data that is supplied by the
user. data
is a pointer to the data to be passed to the user-supplied function. See
the “Introduction, Passing
Data to User-Supplied Functions” at the beginning of this manual for
more details.
Description
Function imsls_f_nonlinear_regression fits a nonlinear regression model using least squares. The nonlinear regression model is
yi = f(xi; θ) + ɛi i = 1, 2, ..., n
where the observed values of the yi's constitute the responses or values of the dependent variable, the known xi's are the vectors of the values of the independent (explanatory) variables, θ is the vector of p regression parameters, and the ɛi's are independently distributed normal errors with mean 0 and variance σ2. For this model, a least-squares estimate of θ is also a maximum likelihood estimate of θ.
The residuals for the model are as follows:
ei(θ) = yi – f(xi; θ) i = 1, 2, ..., n
A value of θ that minimizes
is a least-squares estimate of θ. Function imsls_f_nonlinear_regression is designed so that the values of the function f(xi; θ) are computed one at a time by a user-supplied function.
Function imsls_f_nonlinear_regression is based on MINPACK routines LMDIF and LMDER by Moré et al. (1980) that use a modified Levenberg-Marquardt method to generate a sequence of approximations to a minimum point. Let
be the current estimate of θ. A new estimate is given by
where sc is a solution to the following:
Here
is the Jacobian evaluated at
The algorithm uses a “trust region” approach with a step bound of δc. A solution of the equations is first obtained for
μc = 0. If ||sc||2 < δc
this update is accepted; otherwise, μc is set to a positive value and another solution is obtained. The method is discussed by Levenberg (1944), Marquardt (1963), and Dennis and Schnabel (1983, pp. 129−147, 218−338).
If a user-supplied function is specified in IMSLS_JACOBIAN, the Jacobian is computed analytically; otherwise, forward finite differences are used to estimate the Jacobian numerically. In the latter case, especially if type float is used, the estimate of the Jacobian may be so poor that the algorithm terminates at a noncritical point. In such instances, the user should either supply a Jacobian function, use type double, or do both.
On some platforms, imsls_f_nonlinear_regression can evaluate the user-supplied functions fcn and jacobian in parallel. This is done only if the function imsls_omp_options is called to flag user-defined functions as thread-safe. A function is thread-safe if there are no dependencies between calls. Such dependencies are usually the result of writing to global or static variables.
Programming Notes
Nonlinear regression allows substantial flexibility over linear regression because the user can specify the functional form of the model. This added flexibility can cause unexpected convergence problems for users that are unaware of the limitations of the software. Also, in many cases, there are possible remedies that may not be immediately obvious. The following is a list of possible convergence problems and some remedies. There is not a one-to-one correspondence between the problems and the remedies. Remedies for some problems also may be relevant for the other problems.
1. A local minimum is found. Try a different starting value. Good starting values often can be obtained by fitting simpler models. For example, for a nonlinear function
good starting values can be obtained from the estimated linear regression coefficients
and
from a simple linear regression of ln y on ln x. The starting values for the nonlinear regression in this case would be
If an approximate linear model is not clear, then simplify the model by reducing the number of nonlinear regression parameters. For example, some nonlinear parameters for which good starting values are known could be set to these values in order to simplify the model for computing starting values for the remaining parameters.
2. The estimate of θ is incorrectly returned as the same or very close to the initial estimate. This occurs often because of poor scaling of the problem, which might result in the residual sum of squares being either very large or very small relative to the precision of the computer. The optional arguments allow control of the scaling.
3. The model is discontinuous as a function of θ. (The function f(x;θ) can be a discontinuous function of x.)
4. Overflow occurs during the computations. Make sure the user-supplied functions do not overflow at some value of θ.
5. The estimate of θ is going to infinity. A parameterization of the problem in terms of reciprocals may help.
6. Some components of θ are outside known bounds. This can sometimes be handled by making a function that produces artificially large residuals outside of the bounds (even though this introduces a discontinuity in the model function).
Examples
Example 1
In this example (Draper and Smith 1981, p. 518), the following nonlinear model is fit:
#include <math.h>
#include <imsls.h>
float fcn(int, float[], int, float[]);
int main ()
{
#define N_OBSERVATIONS 4
int n_independent = 1;
int n_parameters = 2;
float *theta_hat;
float x[N_OBSERVATIONS][1] = {10.0, 20.0, 30.0, 40.0};
float y[N_OBSERVATIONS] = {0.48, 0.42, 0.40, 0.39};
imsls_omp_options(IMSLS_SET_FUNCTIONS_THREAD_SAFE, 1, 0);
/* Nonlinear regression */
theta_hat = imsls_f_nonlinear_regression(fcn, n_parameters,
N_OBSERVATIONS, n_independent, (float *)x, y, 0);
/* Print estimates */
imsls_f_write_matrix("estimated coefficients", 1, n_parameters,
theta_hat, 0);
} /* End of main */
float fcn(int n_independent, float x[], int n_parameters, float theta[])
{
return (theta[0] + (0.49 - theta[0])*exp(theta[1]*(x[0] - 8)));
} /* End of fcn */
Output
estimated coefficients
1 2
0.3807 -0.0794
Example 2
Consider the nonlinear regression model and data set discussed by Neter et al. (1983, pp. 475−478):
There are two parameters and one independent variable. The data set considered consists of 15 observations.
#include <imsls.h>
#include <imsl.h>
#include <stdio.h>
#include <math.h>
static float fcn(int, float[], int, float[]);
static void jacobian(int, float[], int, float[], float[]);
int main()
{
int df;
int n_independent=1;
int n_parameters = 2;
int n_obs = 15;
float *theta_hat, *r, *y_hat, *var_covar;
float grad_eps = 1.0e-9;
float theta_guess[2] = {60.0, -0.03};
float a, dfe, normalValue;
float y[15] =
{
54.0, 50.0, 45.0, 37.0, 35.0,
25.0, 20.0, 16.0, 18.0, 13.0,
8.0, 11.0, 8.0, 4.0, 6.0
};
float x[15] =
{
2.0, 5.0, 7.0, 10.0, 14.0,
19.0, 26.0, 31.0, 34.0, 38.0,
45.0, 52.0, 53.0, 60.0, 65.0
};
char *fmt="%15.3f";
char *dashes=
"-------------------------------------------------------------";
/* Nonlinear regression */
theta_hat = imsls_f_nonlinear_regression(fcn, n_parameters,
n_obs, n_independent, x, y,
IMSLS_THETA_GUESS, theta_guess,
IMSLS_GRADIENT_EPS, grad_eps,
IMSLS_DF, &df,
IMSLS_R, &r,
IMSLS_PREDICTED, &y_hat,
IMSLS_VARIANCE_COVARIANCE_MATRIX, &var_covar,
IMSLS_JACOBIAN, jacobian, 0);
/* Print results */
/* Calculate and Print Coefficients & their 95% Confidence Limits */
printf(" \n ESTIMATED COEFFICIENTS \n");
printf("%s\n", dashes);
printf(" Coefficient | Lower 95%% Limit |");
printf(" Estimate | Upper 95%% Limit \n");
dfe = (float) df;
normalValue = imsls_f_t_inverse_cdf(0.975, dfe);
a = normalValue * sqrt(var_covar[0]);
printf(" Theta_1 | %10.3f | %7.3f | %12.3f\n",
theta_hat[0] - a, theta_hat[0], theta_hat[0] + a);
a = normalValue * sqrt(var_covar[3]);
printf(" Theta_2 | %10.3f | %7.3f | %12.3f\n",
theta_hat[1] - a, theta_hat[1], theta_hat[1] + a);
printf("%s\n", dashes);
imsls_f_write_matrix("Var/Covar matrix", n_parameters, n_parameters,
var_covar, IMSLS_WRITE_FORMAT, fmt, 0);
imsls_f_write_matrix("Predicted values", 1, n_obs, y_hat,
IMSLS_WRITE_FORMAT, "%7.2f", 0);
return 0;
}
static float fcn(int n_independent, float x[],
int n_parameters, float theta[])
{
return (theta[0]*exp(x[0]*theta[1]));
} /* End of fcn */
static void jacobian(int n_independent, float x[],
int n_parameters, float theta[],
float fjac[])
{
fjac[0] = exp(theta[1]*x[0]);
fjac[1] = theta[0]*x[0]*exp(theta[1]*x[0]);
} /* End of jacobian */
Output
ESTIMATED COEFFICIENTS
-------------------------------------------------------------
Coefficient | Lower 95% Limit | Estimate | Upper 95% Limit
Theta_1 | 55.426 | 58.607 | 61.787
Theta_2 | -0.043 | -0.040 | -0.036
-------------------------------------------------------------
Var/Covar matrix
1 2
1 2.167 -0.002
2 -0.002 0.000
Predicted values
1 2 3 4 5 6 7 8
54.15 48.08 44.42 39.45 33.67 27.62 20.94 17.18
9 10 11 12 13 14 15
15.26 13.02 9.87 7.48 7.19 5.45 4.40
Informational Errors
IMSLS_STEP_TOLERANCE Scaled step tolerance satisfied. The current point may be an approximate local solution, but it is also possible that the algorithm is making very slow progress and is not near a solution or that “step_eps” is too big.
Warning Errors
IMSLS_LITTLE_FCN_CHANGE Both the actual and predicted relative reductions in the function are less than or equal to the relative function tolerance.
IMSLS_TOO_MANY_ITN Maximum number of iterations exceeded.
IMSLS_TOO_MANY_JACOBIAN_EVAL Maximum number of Jacobian evaluations exceeded.
IMSLS_UNBOUNDED Five consecutive steps have been taken with the maximum step length.
IMSLS_FALSE_CONVERGENCE The iterates appear to be converging to a noncritical point.
Fatal Errors
IMSLS_TOO_MANY_FCN_EVAL Maximum number of function evaluations exceeded.
|
Visual Numerics, Inc. PHONE: 713.784.3131 FAX:713.781.9260 |