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,
0)
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 = "#". |