nonlinear_regression

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,

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 θ.

Convergence is declared if |gi| * max{|θi|, 1.0/ti}/SSE is less than gradient_eps for i= 0, 1, 2, …, n_parameters, where gi is the i‑th component of an internal intermediate gradient vector.
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

Convergence is declared if |gn+ i|/ max{|θi|, 1.0/ti} is less than step_eps for i = 0, 1, 2, …, n, where gn+ i is the i‑th component of the last step and n = n_parameters.

Default: step_eps = ɛ2/3,where ɛ is the machine precision

IMSLS_SSE_REL_EPS, float sse_rel_eps (Input)
Relative SSE function tolerance.

Convergence is declared if the change in SSE is less than or equal to sse_rel_eps * SSE in absolute value.

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.

Convergence is declared if SSE is less than sse_abs_eps.

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 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 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

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:

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.

The first stopping criterion for imsls_f_nonlinear_regression occurs when SSE is less than the absolute function tolerance. The second stopping criterion occurs when the norm of the scaled gradient is less than the given gradient tolerance. The third stopping criterion occurs when the scaled distance between the last two steps is less than the step tolerance. The third stopping criterion also generates error IMSLS_LITTLE_FCN_CHANGE. The fourth stopping criterion occurs when the relative change in SSE is less than sse_rel_eps. The fourth stopping criterion also generates error code IMSLS_FALSE_CONVERGENCE. See Dennis and Schnabel (1983, pages 159−161, 278−280, and 347−348) for a discussion of stopping criteria and choice of tolerances.

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 functionf(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);

}

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)));

}

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 <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);

}

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.
IMSLS_STOP_USER_FCN	Request from user supplied function to stop algorithm. User flag = "#".