NCL Home > Documentation > Functions > General applied math

# lspoly

Calculates a set of coefficients for a weighted least squares polynomial fit to the given data.

## Prototype

function lspoly (
x     : numeric,
y     : numeric,
wgt   : numeric,
n [1] : integer
)

return_val  :  float or double

## Arguments

x

Abscissa values of the data.

This can be one-dimensional or multi-dimensional. If one-dimensional, it must be the same length as the rightmost dimension of y. If multi-dimensional, it must be the same dimensionality as y.

y

Ordinate values of the data.

This can be one-dimensional or multi-dimensional. If one-dimensional, it must be the same length as the rightmost dimension of x. If multi-dimensional, it must be the same dimensionality as x.

wgt

Weights for a weighted least squares model. If all data values are to be assigned equal weights, then setting the argument equal to a scalar 1.0 will result in all the weights being set to 1.0. Note: if x or y is equal to _FillValue (if present), the weight will be set to 0.0 for that coordinate pair.

n

The number of coefficients desired (i.e., n-1 will be the degree of the polynomial). Due to the computational method used, n should be less than or equal to five.

## Return value

The return array will have the same dimensionality as y, except the rightmost dimension will be of length n.

If either x or y are of type double, then the return array is returned as double. Otherwise, the returned coefficients are returned as type float.

## Description

Given a set of data (x(i),y(i)), i = 1,...,m, lspoly calculates a set of coefficients for a weighted least squares polynomial fit to the given data. It is necessary that the number of data points) be greater than n (the number of coefficients).

Accuracy: for lower order polynomials (n .le. 5), lspoly can be expected to give satisfactory results.

Algorithm: lspoly forms the normal and solves the resulting square linear system using gaussian elimination with full pivoting.

Note: In NCL versions 6.1.2 and earlier, this function would fail in certain cases where the input data were "scaled" down. We replaced this routine with a SLATEC version. It takes the same arguments, but uses a different algorithm under the hood. You can get the "old" (pre NCL V6.2.0) version of lspoly by using "lspoly_old", but we don't recommend it for regular use.

Use lspoly_n if you want to specify which dimension(s) to do the calculation across.

## Examples

Example 1

x = (/-4.5, -3.2, -1.4, 0.8, 2.5, 4.1/)
y = (/ 0.7,  2.3,  3.8, 5.0, 5.5, 5.6/)

n = 4
c = lspoly(x,y, 1, n)    ; all weights are set to one
print(c)
The 3rd degree polynomial is
Y = c(0) + c(1)*x + c(2)*x^2 + c(3)*x^3
The coefficients (which agree with those returned from Mathematica) are:
(0)      4.66863
(1)      0.489392
(2)     -0.0742387
(3)      0.00267663

Example 3: Compare with an example from Matlab's polyfit: Fit Polynomial to Error Function.

diri = "./"
fili = "polyfit.matlab.txt"
pthi = diri+fili

nrow = 16
ncol = 4
x    = data(:,0)
y    = data(:,1)
matfit = data(:,2)             ; polyfit
materr = data(:,3)             ; diff: (y-matfit)

n = 7
c = lspoly(x, y, 1, n)         ; wgt=1
print(c)

nclfit = c(0) + c(1)*x   + c(2)*x^2 + c(3)*x^3 \
+ c(4)*x^4 + c(5)*x^5 + c(6)*x^6
nclerr = y - nclfit

print(sprintf("%5.1f", x)+"  "+sprintf("%9.5f", y)   \
+"  "+sprintf("%9.5f", matfit)+"  "+sprintf("%9.5f", materr) \
+"  "+sprintf("%9.5f", nclfit)+"  "+sprintf("%9.5f", nclerr))

msqerr = avg(materr^2)
nsqerr = avg(nclerr^2)

print(sprintf("%12.10f", msqerr)+"  "+sprintf("%12.10f", nsqerr) )

The (edited) output is:

Variable: c
Type: float

Dimensions and sizes:        [7]        <== n=7

(0)  1.240614e-05
(1)  1.12672
(2)  0.01681045
(3)  -0.435255
(4)  0.08767352
(5)  0.07092082
(6)  -0.02416329

Matlab     Matlab       NCL       NCL
X        Y       polyfit     diff       lspoly     diff

(0)    0.0    0.00000    0.00044   -0.00044    0.00001   -0.00001
(1)    0.1    0.11246    0.11185    0.00061    0.11243    0.00003
(2)    0.2    0.22270    0.22231    0.00039    0.22271   -0.00001
(3)    0.3    0.32863    0.32872   -0.00010    0.32865   -0.00002
(4)    0.4    0.42839    0.42880   -0.00041    0.42841   -0.00002
(5)    0.5    0.52050    0.52093   -0.00043    0.52049    0.00001
(6)    0.6    0.60386    0.60408   -0.00023    0.60383    0.00003
(7)    0.7    0.67780    0.67775    0.00005    0.67779    0.00001
(8)    0.8    0.74210    0.74183    0.00027    0.74211   -0.00001
(9)    0.9    0.79691    0.79654    0.00037    0.79693   -0.00002
(10)   1.0    0.84270    0.84238    0.00032    0.84272   -0.00002
(11)   1.1    0.88021    0.88005    0.00016    0.88020    0.00001
(12)   1.2    0.91031    0.91035   -0.00004    0.91028    0.00003
(13)   1.3    0.93401    0.93422   -0.00021    0.93400    0.00001
(14)   1.4    0.95229    0.95258   -0.00030    0.95233   -0.00004
(15)   1.5    0.96611    0.96639   -0.00028    0.96610    0.00001

Mean Square Fit Difference
polyfit       lspoly

(0)  0.0000001047  0.0000000004