hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_fit_1dcheb_arb (e02ad)

Purpose

nag_fit_1dcheb_arb (e02ad) computes weighted least squares polynomial approximations to an arbitrary set of data points.

Syntax

[a, s, ifail] = e02ad(kplus1, x, y, w, 'm', m)
[a, s, ifail] = nag_fit_1dcheb_arb(kplus1, x, y, w, 'm', m)

Description

nag_fit_1dcheb_arb (e02ad) determines least squares polynomial approximations of degrees 0,1,,k0,1,,k to the set of data points (xr,yr)(xr,yr) with weights wrwr, for r = 1,2,,mr=1,2,,m.
The approximation of degree ii has the property that it minimizes σiσi the sum of squares of the weighted residuals εrεr, where
εr = wr(yrfr)
εr=wr(yr-fr)
and frfr is the value of the polynomial of degree ii at the rrth data point.
Each polynomial is represented in Chebyshev series form with normalized argument xx-. This argument lies in the range 1-1 to + 1+1 and is related to the original variable xx by the linear transformation
x = ((2xxmaxxmin))/((xmaxxmin)).
x-= (2x-xmax-xmin) (xmax-xmin) .
Here xmaxxmax and xminxmin are respectively the largest and smallest values of xrxr. The polynomial approximation of degree ii is represented as
(1/2)ai + 1,1T0(x) + ai + 1,2T1(x) + ai + 1,3T2(x) + + ai + 1,i + 1Ti(x),
12ai+1,1T0(x-)+ai+1,2T1(x-)+ai+1,3T2(x-)++ai+1,i+1Ti(x-),
where Tj(x)Tj(x-), for j = 0,1,,ij=0,1,,i, are the Chebyshev polynomials of the first kind of degree jj with argument (x)(x-).
For i = 0,1,,ki=0,1,,k, the function produces the values of ai + 1,j + 1ai+1,j+1, for j = 0,1,,ij=0,1,,i, together with the value of the root-mean-square residual si = sqrt(σi / (mi1))si=σi/(m-i-1). In the case m = i + 1m=i+1 the function sets the value of sisi to zero.
The method employed is due to Forsythe (1957) and is based on the generation of a set of polynomials orthogonal with respect to summation over the normalized dataset. The extensions due to Clenshaw (1960) to represent these polynomials as well as the approximating polynomials in their Chebyshev series forms are incorporated. The modifications suggested by Reinsch and Gentleman (see Gentleman (1969)) to the method originally employed by Clenshaw for evaluating the orthogonal polynomials from their Chebyshev series representations are used to give greater numerical stability.
For further details of the algorithm and its use see Cox (1974) and Cox and Hayes (1973).
Subsequent evaluation of the Chebyshev series representations of the polynomial approximations should be carried out using nag_fit_1dcheb_eval (e02ae).

References

Clenshaw C W (1960) Curve fitting with a digital computer Comput. J. 2 170–173
Cox M G (1974) A data-fitting package for the non-specialist user Software for Numerical Mathematics (ed D J Evans) Academic Press
Cox M G and Hayes J G (1973) Curve fitting: a guide and suite of algorithms for the non-specialist user NPL Report NAC26 National Physical Laboratory
Forsythe G E (1957) Generation and use of orthogonal polynomials for data fitting with a digital computer J. Soc. Indust. Appl. Math. 5 74–88
Gentleman W M (1969) An error analysis of Goertzel's (Watt's) method for computing Fourier coefficients Comput. J. 12 160–165
Hayes J G (ed.) (1970) Numerical Approximation to Functions and Data Athlone Press, London

Parameters

Compulsory Input Parameters

1:     kplus1 – int64int32nag_int scalar
k + 1k+1, where kk is the maximum degree required.
Constraint: 0 < kplus1mdist0<kplus1mdist, where mdistmdist is the number of distinct xx values in the data.
2:     x(m) – double array
m, the dimension of the array, must satisfy the constraint mmdist2mmdist2, where mdistmdist is the number of distinct xx values in the data.
The values xrxr of the independent variable, for r = 1,2,,mr=1,2,,m.
Constraint: the values must be supplied in nondecreasing order with x(m) > x(1)xm>x1.
3:     y(m) – double array
m, the dimension of the array, must satisfy the constraint mmdist2mmdist2, where mdistmdist is the number of distinct xx values in the data.
The values yryr of the dependent variable, for r = 1,2,,mr=1,2,,m.
4:     w(m) – double array
m, the dimension of the array, must satisfy the constraint mmdist2mmdist2, where mdistmdist is the number of distinct xx values in the data.
The set of weights, wrwr, for r = 1,2,,mr=1,2,,m. For advice on the choice of weights, see Section [Weighting of data points] in the E02 Chapter Introduction.
Constraint: w(r) > 0.0wr>0.0, for r = 1,2,,mr=1,2,,m.

Optional Input Parameters

1:     m – int64int32nag_int scalar
Default: The dimension of the arrays x, y, w. (An error is raised if these dimensions are not equal.)
The number mm of data points.
Constraint: mmdist2mmdist2, where mdistmdist is the number of distinct xx values in the data.

Input Parameters Omitted from the MATLAB Interface

lda work1 work2

Output Parameters

1:     a(lda,kplus1) – double array
ldakplus1ldakplus1.
The coefficients of Tj(x)Tj(x-) in the approximating polynomial of degree ii. a(i + 1,j + 1)ai+1j+1 contains the coefficient ai + 1,j + 1ai+1,j+1, for i = 0,1,,ki=0,1,,k and j = 0,1,,ij=0,1,,i.
2:     s(kplus1) – double array
s(i + 1)si+1 contains the root-mean-square residual sisi, for i = 0,1,,ki=0,1,,k, as described in Section [Description]. For the interpretation of the values of the sisi and their use in selecting an appropriate degree, see Section [General] in the E02 Chapter Introduction.
3:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
The weights are not all strictly positive.
  ifail = 2ifail=2
The values of x(r)xr, for r = 1,2,,mr=1,2,,m, are not in nondecreasing order.
  ifail = 3ifail=3
All x(r)xr have the same value: thus the normalization of x is not possible.
  ifail = 4ifail=4
On entry,kplus1 < 1kplus1<1 (so the maximum degree required is negative)
orkplus1 > mdistkplus1>mdist, where mdistmdist is the number of distinct xx values in the data (so there cannot be a unique solution for degree k = kplus11k=kplus1-1).
  ifail = 5ifail=5
lda < kplus1lda<kplus1.

Accuracy

No error analysis for the method has been published. Practical experience with the method, however, is generally extremely satisfactory.

Further Comments

The time taken is approximately proportional to m(k + 1)(k + 11)m(k+1)(k+11).
The approximating polynomials may exhibit undesirable oscillations (particularly near the ends of the range) if the maximum degree kk exceeds a critical value which depends on the number of data points mm and their relative positions. As a rough guide, for equally-spaced data, this critical value is about 2 × sqrt(m)2×m. For further details see page 60 of Hayes (1970).

Example

function nag_fit_1dcheb_arb_example
kplus1 = int64(4);
x = [1;
     2.1;
     3.1;
     3.9;
     4.9;
     5.8;
     6.5;
     7.1;
     7.8;
     8.4;
     9];
y = [10.4;
     7.9;
     4.7;
     2.5;
     1.2;
     2.2;
     5.1;
     9.2;
     16.1;
     24.5;
     35.3];
w = [1;
     1;
     1;
     1;
     1;
     0.8;
     0.8;
     0.7;
     0.5;
     0.3;
     0.2];
[a, s, ifail] = nag_fit_1dcheb_arb(kplus1, x, y, w)
 

a =

   12.1740         0         0         0
   12.2954    0.2740         0         0
   20.7345    6.2016    8.1876         0
   24.1429    9.4065   10.8400    3.0589


s =

    4.0659
    4.2840
    1.6865
    0.0682


ifail =

                    0


function e02ad_example
kplus1 = int64(4);
x = [1;
     2.1;
     3.1;
     3.9;
     4.9;
     5.8;
     6.5;
     7.1;
     7.8;
     8.4;
     9];
y = [10.4;
     7.9;
     4.7;
     2.5;
     1.2;
     2.2;
     5.1;
     9.2;
     16.1;
     24.5;
     35.3];
w = [1;
     1;
     1;
     1;
     1;
     0.8;
     0.8;
     0.7;
     0.5;
     0.3;
     0.2];
[a, s, ifail] = e02ad(kplus1, x, y, w)
 

a =

   12.1740         0         0         0
   12.2954    0.2740         0         0
   20.7345    6.2016    8.1876         0
   24.1429    9.4065   10.8400    3.0589


s =

    4.0659
    4.2840
    1.6865
    0.0682


ifail =

                    0



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013