naginterfaces.library.fit.dim2_​spline_​panel

naginterfaces.library.fit.dim2_spline_panel(x, y, f, w, lamda, mu, point, eps)

dim2_spline_panel forms a minimal, weighted least squares bicubic spline surface fit with prescribed knots to a given set of data points.

For full information please refer to the NAG Library document for e02da

https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02daf.html

Parameters

xfloat, array-like, shape $\left(m\right)$

The coordinates of the data point $(x_{\text{r}},y_{\text{r}},f_{\text{r}})$, for $\text{r}=1,2,\dots ,m$. The order of the data points is immaterial, but see the array $point$.

yfloat, array-like, shape $\left(m\right)$

The coordinates of the data point $(x_{\text{r}},y_{\text{r}},f_{\text{r}})$, for $\text{r}=1,2,\dots ,m$. The order of the data points is immaterial, but see the array $point$.

ffloat, array-like, shape $\left(m\right)$

The coordinates of the data point $(x_{\text{r}},y_{\text{r}},f_{\text{r}})$, for $\text{r}=1,2,\dots ,m$. The order of the data points is immaterial, but see the array $point$.

wfloat, array-like, shape $\left(m\right)$

The weight $w_r$ of the $r$th data point. It is important to note the definition of weight implied by the equation (1) in Notes, since it is also common usage to define weight as the square of this weight. In this function, each $w_r$ should be chosen inversely proportional to the (absolute) accuracy of the corresponding $f_r$, as expressed, for example, by the standard deviation or probable error of the $f_r$. When the $f_r$ are all of the same accuracy, all the $w_r$ may be set equal to $1.0$.

lamdafloat, array-like, shape $\left(\text{px}\right)$

$lamda[\text{i}+3]$ must contain the $\text{i}$th interior knot ${\lambda }_{\text{i}+4}$ associated with the variable $x$, for $\text{i}=1,2,\dots ,\text{px}-8$. The knots must be in nondecreasing order and lie strictly within the range covered by the data values of $x$. A knot is a value of $x$ at which the spline is allowed to be discontinuous in the third derivative with respect to $x$, though continuous up to the second derivative. This degree of continuity can be reduced, if you require, by the use of coincident knots, provided that no more than four knots are chosen to coincide at any point. Two, or three, coincident knots allow loss of continuity in, respectively, the second and first derivative with respect to $x$ at the value of $x$ at which they coincide. Four coincident knots split the spline surface into two independent parts. For choice of knots see Further Comments.

mufloat, array-like, shape $\left(\text{py}\right)$

$mu[\text{i}+3]$ must contain the $\text{i}$th interior knot ${\mu }_{\text{i}+4}$ associated with the variable $y$, for $\text{i}=1,2,\dots ,\text{py}-8$.

pointint, array-like, shape $\left(\text{npoint}\right)$

Indexing information usually provided by dim2_spline_sort() which enables the data points to be accessed in the order which produces the advantageous matrix structure mentioned in Notes. This order is such that, if the $(x,y)$ plane is thought of as being divided into rectangular panels by the two sets of knots, all data in a panel occur before data in succeeding panels, where the panels are numbered from bottom to top and then left to right with the usual arrangement of axes, as indicated in Figure [label omitted].

[figure omitted]

A data point lying exactly on one or more panel sides is considered to be in the highest numbered panel adjacent to the point. dim2_spline_sort() should be called to obtain the array $point$, unless it is provided by other means.

epsfloat

A threshold $ϵ$ for determining the effective rank of the system of linear equations. The rank is determined as the number of elements of the array $dl$ which are nonzero. An element of $dl$ is regarded as zero if it is less than $ϵ$. Machine precision is a suitable value for $ϵ$ in most practical applications which have only $2$ or $3$ decimals accurate in data. If some coefficients of the fit prove to be very large compared with the data ordinates, this suggests that $ϵ$ should be increased so as to decrease the rank. The array $dl$ will give a guide to appropriate values of $ϵ$ to achieve this, as well as to the choice of $ϵ$ in other cases where some experimentation may be needed to determine a value which leads to a satisfactory fit.

Returns

lamdafloat, ndarray, shape $\left(\text{px}\right)$

The interior knots $lamda\left[4\right]$ to $lamda[\text{px}-5]$ are unchanged, and the segments $lamda[0:4]$ and $lamda[\text{px}-4:\text{px}]$ contain additional (exterior) knots introduced by the function in order to define the full set of B-splines required. The four knots in the first segment are all set equal to the lowest data value of $x$ and the other four additional knots are all set equal to the highest value: there is experimental evidence that coincident end-knots are best for numerical accuracy. The complete array must be left undisturbed if dim2_spline_evalv() or dim2_spline_evalm() is to be used subsequently.

mufloat, ndarray, shape $\left(\text{py}\right)$

The same remarks apply to $mu$ as to $lamda$ above, with $y$ replacing $x$, and $y$ replacing $x$.

dlfloat, ndarray, shape $((\text{px}-4)\times (\text{py}-4))$

Gives the squares of the diagonal elements of the reduced triangular matrix, divided by the mean squared weight. It includes those elements, less than $ϵ$, which are treated as zero (see Notes).

cfloat, ndarray, shape $((\text{px}-4)\times (\text{py}-4))$

Gives the coefficients of the fit. $c[(\text{py}-4)\times (\text{i}-1)+\text{j}-1]$ is the coefficient $c_{\text{i}\text{j}}$ of Notes and Further Comments, for $\text{j}=1,2,\dots ,\text{py}-4$, for $\text{i}=1,2,\dots ,\text{px}-4$. These coefficients are used by dim2_spline_evalv() or dim2_spline_evalm() to calculate values of the fitted function.

sigmafloat

$\Sigma$, the weighted sum of squares of residuals. This is not computed from the individual residuals but from the right-hand sides of the orthogonally-transformed linear equations. For further details see page 97 of Hayes and Halliday (1974). The two methods of computation are theoretically equivalent, but the results may differ because of rounding error.

rankint

The rank of the system as determined by the value of the threshold $ϵ$.

$rank=(\text{px}-4)\times (\text{py}-4)$

The least squares solution is unique.

$rank\ne (\text{px}-4)\times (\text{py}-4)$

The minimal least squares solution is computed.

Raises

NagValueError

(errno $1$)

At least one set of knots is not in nondecreasing order.

(errno $2$)

More than four knots coincide at a single point.

(errno $3$)

Array $point$ does not indicate the data points in panel order.

(errno $4$)

On entry, $\text{npoint}=⟨value⟩$, $m=⟨value⟩$, $\text{px}=⟨value⟩$ and $\text{py}=⟨value⟩$.

Constraint: $\text{npoint}\ge m+(\text{px}-7)\times (\text{py}-7)$.

(errno $4$)

On entry, $\text{nc}=⟨value⟩$, $\text{px}=⟨value⟩$ and $\text{py}=⟨value⟩$.

Constraint: $\text{nc}=(\text{px}-4)\times (\text{py}-4)$.

(errno $4$)

On entry, $\text{py}=⟨value⟩$.

Constraint: $\text{py}\ge 8$.

(errno $4$)

On entry, $\text{px}=⟨value⟩$.

Constraint: $\text{px}\ge 8$.

(errno $4$)

On entry, $m=⟨value⟩$.

Constraint: $m>1$.

(errno $5$)

All the weights are zero, or rank determined as zero.

Notes

dim2_spline_panel determines a bicubic spline fit $s(x,y)$ to the set of data points $(x_r,y_r,f_r)$ with weights $w_r$, for $\text{r}=1,2,\dots ,m$. The two sets of internal knots of the spline, $\left\{\lambda \right\}$ and $\left\{\mu \right\}$, associated with the variables $x$ and $y$ respectively, are prescribed by you. These knots can be thought of as dividing the data region of the $(x,y)$ plane into panels (see Figure [label omitted] in Parameters). A bicubic spline consists of a separate bicubic polynomial in each panel, the polynomials joining together with continuity up to the second derivative across the panel boundaries.

$s(x,y)$ has the property that $\Sigma$, the sum of squares of its weighted residuals ${\rho }_r$, for $\text{r}=1,2,\dots ,m$, where

$${\rho }_r=w_r(s(x_r,y_r)-f_r)$$

is as small as possible for a bicubic spline with the given knot sets. The function produces this minimized value of $\Sigma$ and the coefficients $c_{ij}$ in the B-spline representation of $s(x,y)$ – see Further Comments. dim2_spline_evalv(), dim2_spline_evalm() and dim2_spline_derivm() are available to compute values and derivatives of the fitted spline from the coefficients $c_{ij}$.

The least squares criterion is not always sufficient to determine the bicubic spline uniquely: there may be a whole family of splines which have the same minimum sum of squares. In these cases, the function selects from this family the spline for which the sum of squares of the coefficients $c_{ij}$ is smallest: in other words, the minimal least squares solution. This choice, although arbitrary, reduces the risk of unwanted fluctuations in the spline fit. The method employed involves forming a system of $m$ linear equations in the coefficients $c_{ij}$ and then computing its least squares solution, which will be the minimal least squares solution when appropriate. The basis of the method is described in Hayes and Halliday (1974). The matrix of the equation is formed using a recurrence relation for B-splines which is numerically stable (see Cox (1972) and de Boor (1972) – the former contains the more elementary derivation but, unlike de Boor (1972), does not cover the case of coincident knots). The least squares solution is also obtained in a stable manner by using orthogonal transformations, viz. a variant of Givens rotation (see Gentleman (1973)). This requires only one row of the matrix to be stored at a time. Advantage is taken of the stepped-band structure which the matrix possesses when the data points are suitably ordered, there being at most sixteen nonzero elements in any row because of the definition of B-splines. First the matrix is reduced to upper triangular form and then the diagonal elements of this triangle are examined in turn. When an element is encountered whose square, divided by the mean squared weight, is less than a threshold $ϵ$, it is replaced by zero and the rest of the elements in its row are reduced to zero by rotations with the remaining rows. The rank of the system is taken to be the number of nonzero diagonal elements in the final triangle, and the nonzero rows of this triangle are used to compute the minimal least squares solution. If all the diagonal elements are nonzero, the rank is equal to the number of coefficients $c_{ij}$ and the solution obtained is the ordinary least squares solution, which is unique in this case.

References

Cox, M G, 1972, The numerical evaluation of B-splines, J. Inst. Math. Appl. (10), 134–149

de Boor, C, 1972, On calculating with B-splines, J. Approx. Theory (6), 50–62

Gentleman, W M, 1973, Least squares computations by Givens transformations without square roots, J. Inst. Math. Applic. (12), 329–336

Hayes, J G and Halliday, J, 1974, The least squares fitting of cubic spline surfaces to general data sets, J. Inst. Math. Appl. (14), 89–103