user_fcn_least

user_fcn_least_squares

Computes a least-squares fit using user-supplied functions.

Synopsis

#include <imsl.h>

float *imsl_f_user_fcn_least_squares (float fcn (int k, float x), int nbasis, int ndata, float xdata[], float ydata[], …, 0)

The type double function is imsl_d_user_fcn_least_squares.

Required Arguments

float fcn (int k, float x) (Input)
User-supplied function that defines the subspace from which the least-squares fit is to be performed. The k-th basis function evaluated at x is f(k, x) where k = 1, 2, …, nbasis.

int nbasis (Input)
Number of basis functions.

int ndata (Input)
Number of data points.

float xdata[] (Input)
Array with ndata components containing the abscissas of the least-squares problem.

float ydata[] (Input)
Array with ndata components containing the ordinates of the least-squares problem.

Return Value

A pointer to the vector containing the coefficients of the basis functions. If a fit cannot be computed, then NULL is returned. To release this space, use imsl_free.

Synopsis with Optional Arguments

#include <imsl.h>

float *imsl_f_user_fcn_least_squares (float fcn(), int nbasis, int ndata, float xdata[], float ydata[],

IMSL_RETURN_USER, float coef[],

IMSL_INTERCEPT, float *intercept,

IMSL_SSE, float *ssq_err,

IMSL_WEIGHTS, float weights[],

IMSL_FCN_W_DATA, float fcn(), void *data,

Optional Arguments

IMSL_RETURN_USER, float coef[] (Output)
The coefficients are stored in the user-supplied array.

IMSL_INTERCEPT, float *intercept (Output)
This option adds an intercept to the model. Thus, the least-squares fit is computed using the user-supplied basis functions augmented by the constant function. The coefficient of the constant function is stored in intercept.

IMSL_SSE, float *ssq_err (Output)
This option returns the error sum of squares.

IMSL_WEIGHTS, float weights[] (Input)
This option requires the user to provide the weights.
Default: all weights equal one

IMSL_FCN_W_DATA, float fcn (int k, float x, float *data), void *data, (Input)
User supplied function that defines the subspace from which the least-squares fit is to be performed, 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 Passing Data to User-Supplied Functions in the introduction to this manual for more details.

Description

The function imsl_f_user_fcn_least_squares computes a best least-squares approximation to given univariate data of the form

by M basis functions

(where M = nbasis). In particular, the default for this function returns the coefficients a which minimize

where w = weights, n = ndata, x = xdata, and f = ydata.

If the optional argument IMSL_INTCERCEPT is chosen, then an intercept is placed in the model, and the coefficients a, returned by imsl_f_user_fcn_least_squares, minimize the error sum of squares as indicated below.

Remarks

Note that although the system is linear, for very large problems the Chapter 8 function, imsl_f_nonlin_least_squares, might be better suited. This is because imsl_f_user_fcn_least_squares will gather the matrix entries one at a time by calls to the user‑supplied function. By using the nonlinear solver and supplying the Jacobian the user can sidestep this issue and likely achieve accurate results since the nonlinear solver utilizes regularization and iterative refinement. Example 3 below demonstrates the use of the nonlinear solver imsl_f_nonlin_least_squares as an alternative to imsl_f_user_fcn_least_squares.

Examples

Example 1

This example fits the following two functions (indexed by δ):

1 + sinx + 7 sin3x + δε

where ɛ is a random uniform deviate over the range [‑1, 1] and δ is 0 for the first function and 1 for the second. These functions are evaluated at 90 equally spaced points on the interval [0, 6]. Four basis functions are used: 1, sinx, sin2x, sin3x.

#include <imsl.h>

#include <stdio.h>

#include <math.h>

#define NDATA 90

/* Define function */

#define F(x) (float)(1.+ sin(x)+7.*sin(3.0*x))

float fcn(int n, float x);

int main()

{

int nbasis = 4, i, delta;

float ydata[NDATA], xdata[NDATA], *random, *coef;

/* Generate random numbers */

imsl_random_seed_set(1234567);

random = imsl_f_random_uniform(NDATA, 0);

/* Set up data */

for(delta = 0; delta < 2; delta++) {

for (i = 0; i < NDATA; i++) {

xdata[i] = 6.*(float)i /((float)(NDATA-1));

ydata[i] = F(xdata[i]) + (delta)*2.*(random[i]-.5);

}

coef = imsl_f_user_fcn_least_squares(fcn, nbasis, NDATA, xdata,

ydata, 0);

printf("\nFor delta = %1d", delta);

imsl_f_write_matrix("the computed coefficients are\n",

1, nbasis, coef, 0);

}

float fcn(int n, float x)

{

return (n == 1) ? 1.0 : sin((n-1)*x);

}

Output

For delta = 0

the computed coefficients are

1 2 3 4

1 1 -0 7

For delta = 1

the computed coefficients are

1 2 3 4

0.979 0.998 0.096 6.839

Example 2

Recall that the first example fitted the following two functions (indexed by δ):

1 + sinx + 7 sin3x + δε

where ε is a random uniform deviate over the range[−1, 1] , and δ is 0 for the first function and 1 for the second. These functions are evaluated at 90 equally spaced points on the interval [0, 6]. Previously, the four basis functions were used: 1, sinx, sin2x, sin3x. This example uses the four basis functions: sinx, sin2x, sin3x, sin4x, combined with the intercept option.

#include <imsl.h>

#include <stdio.h>

#include <math.h>

#define NDATA 90

/* Define function */

#define F(x) (float)(1.+ sin(x)+7.*sin(3.0*x))

float fcn(int n, float x);

int main()

{

int nbasis = 4, i, delta;

float ydata[NDATA], xdata[NDATA], *random, *coef, intercept;

/* Generate random numbers */

imsl_random_seed_set(1234567);

random = imsl_f_random_uniform(NDATA, 0);

/* Set up data */

for(delta = 0; delta < 2; delta++){

for (i = 0; i < NDATA; i++) {

xdata[i] = 6.*(float)i /((float)(NDATA-1));

ydata[i] = F(xdata[i]) + (delta)*2.*(random[i]-.5);

}

coef = imsl_f_user_fcn_least_squares(fcn, nbasis, NDATA, xdata,

ydata,

IMSL_INTERCEPT, &intercept,

0);

printf("\nFor delta = %1d\n", delta);

printf("The predicted intercept value is %10.3f\n" ,

intercept);

imsl_f_write_matrix("the computed coefficients are\n",

1, nbasis, coef, 0);

}

float fcn(int n, float x)

{

return sin(n*x);

}

Output

For delta = 0

The predicted intercept value is 1.000

the computed coefficients are

1 2 3 4

1 0 7 -0

For delta = 1

The predicted intercept value is 0.978

the computed coefficients are

1 2 3 4

0.998 0.097 6.841 0.075

Example 3

This example solves the same problem as Example 1, demonstrating the use of imsl_f_nonlin_least_squares as an alternative to imsl_f_user_fcn_least_squares. See the Remarks section above for a discussion of why, for some problems, imsl_f_nonlinear_least_squares might be better suited than imsl_f_user_fcn_least_squares.

#include <stdio.h>

#include <imsl.h>

#include <math.h>

float fcn(int n, float x);

void fcn2(int m, int n, float x[], float f[], void* data);

void fcn2jac(int m, int n, float x[], float fjac[], int fjac_col_dim,

void *data);

* The following structure type will be used to pass

* additional data to imsl_f_nonlin_least_squares() using

* the optional argument IMSL_FCN_W_DATA

typedef struct {

float *xdata;

float *ydata;

} Problem_data;

int main()

{

#define NDATA 90

#define F(x) (float)(1.+ sin(x)+7.*sin(3.0*x))

int nbasis = 4, i, delta;

float *random, *coef, *coef2;

float ydata[NDATA], xdata[NDATA];

Problem_data data;

data.xdata = xdata;

data.ydata = ydata;

imsl_random_seed_set(1234567);

random = imsl_f_random_uniform(NDATA, 0);

for(delta = 0; delta < 2; delta++) {

printf("\n\nFor delta = %1d", delta);

for (i = 0; i < NDATA; i++) {

xdata[i] = 6.*(float)i /((float)(NDATA-1));

ydata[i] = F(xdata[i]) + (delta)*2.*(random[i]-.5);

}

coef = imsl_f_user_fcn_least_squares(fcn, nbasis, NDATA, xdata,

ydata, 0);

imsl_f_write_matrix(

" Coefficients from user_fcn_least_squares\n",

1, nbasis, coef, IMSL_NO_COL_LABELS, IMSL_NO_ROW_LABELS,

IMSL_WRITE_FORMAT, "%6.3f", 0);

coef2 = imsl_f_nonlin_least_squares(NULL, NDATA, nbasis,

IMSL_FCN_W_DATA, fcn2, &data,

IMSL_JACOBIAN_W_DATA, fcn2jac, &data, 0);

imsl_f_write_matrix(

" Coefficients from nonlin_least_squares \n",

1, nbasis, coef2, IMSL_NO_COL_LABELS, IMSL_NO_ROW_LABELS,

IMSL_WRITE_FORMAT, "%6.3f", 0);

}

float fcn(int n, float x)

{

return (n == 1) ? 1.0 : sin((n-1)*x);

}

void fcn2(int m, int n, float x[], float f[], void *data)

{

int i;

float *xdata, *ydata;

xdata = ((Problem_data*)data)->xdata;

ydata = ((Problem_data*)data)->ydata;

for (i=0;i<m;i++)

f[i] = x[0] + x[1]*sin(xdata[i]) + x[2]*sin(2.0*xdata[i]) +

x[3]*sin(3.0*xdata[i]) - ydata[i];

}

void fcn2jac(int m, int n, float x[], float fjac[], int fjac_col_dim,

void *data)

{

int i, j;

float *xdata;

xdata = ((Problem_data*)data)->xdata;

for (i=0;i<m;i++)

for (j=0;j<n;j++)

fjac[i*fjac_col_dim+j] = (j==0) ? 1.0 : sin(j*xdata[i]);

}

Output

For delta = 0

Coefficients from user_fcn_least_squares

1.000 1.000 0.000 7.000

Coefficients from nonlin_least_squares

1.000 1.000 0.000 7.000

For delta = 1

Coefficients from user_fcn_least_squares

0.979 0.998 0.096 6.839

Coefficients from nonlin_least_squares

0.979 0.998 0.096 6.839

Warning Errors

IMSL_LINEAR_DEPENDENCE	Linear dependence of the basis functions exists. One or more components of coef are set to zero.
IMSL_LINEAR_DEPENDENCE_CONST	Linear dependence of the constant function and basis functions exists. One or more components of coef are set to zero.

Fatal Errors

IMSL_NEGATIVE_WEIGHTS_2	All weights must be greater than or equal to zero.
IMSL_STOP_USER_FCN	Request from user supplied function to stop algorithm. User flag = "#".