# Created by Octave 3.4.0, Mon Nov 28 01:51:38 2011 UTC <builder@apiary.egro.altlinux.org>
# name: cache
# type: cell
# rows: 3
# columns: 1
# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
gcvspl


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 10652
B-spline data smoothing subroutine.

[yf,wk]=gcvspl(x(n,1),y(n,k), xf(nf,1)=x, wx(n,1)=[],wy(1,k)=[], m=2,md=2,val=1, ider=[0])

uses  GCVSPL.FOR, 1986-05-12
from  http://www.netlib.org/gcv/index.html
for   B-spline data smoothing using generalized cross-validation
      and mean squared prediction or explicit user smoothing
by    H.J. Woltring,  University of Nijmegen,
      Philips Medical Systems, Eindhoven (The Netherlands)

Purpose:
      Natural B-spline data smoothing subroutine, using the Generali-
      zed Cross-Validation and Mean-Squared Prediction Error Criteria
      of Craven & Wahba (1979). Alternatively, the amount of smoothing
      can be given explicitly, or it can be based on the effective
      number of degrees of freedom in the smoothing process as defined
      by Wahba (1980). The model assumes uncorrelated, additive noise
      and essentially smooth, underlying functions. The noise may be
      non-stationary, and the independent co-ordinates may be spaced
      non-equidistantly. Multiple datasets, with common independent
      variables and weight factors are accomodated.
      A full description of the package is provided in:
      H.J. Woltring (1986), A FORTRAN package for generalized,
      cross-validatory spline smoothing and differentiation.
      Advances in Engineering Software 8(2):104-113

Meaning of parameters:
      X(N,1)  Independent variables: strictly increasing knot
              sequence, with X(I-1).lt.X(I), I=2,...,N.
      Y(N,K)  Input data to be smoothed (or interpolated).
      XF(NF,1)Points where the function should be approximated.
      WX(N,1) Weight factor array; WX(I) corresponds with
              the relative inverse variance of point Y(I,*).
              If no relative weighting information is
              available, the WX(I) should be set to ONE.
              All WX(I).gt.ZERO, I=1,...,N.
      WY(1,K) Weight factor array; WY(J) corresponds with
              the relative inverse variance of point Y(*,J).
              If no relative weighting information is
              available, the WY(J) should be set to ONE.
              All WY(J).gt.ZERO, J=1,...,K.
              NB: The effective weight for point Y(I,J) is
              equal to WX(I)*WY(J).
      M       Half order of the required B-splines (spline
              degree 2*M-1), with M.gt.0. The values M =
              1,2,3,4 correspond to linear, cubic, quintic,
              and heptic splines, respectively. N.ge.2*M.
      MD      Optimization mode switch:
              |MD| = 1: Prior given value for p in VAL
                        (VAL.ge.ZERO). This is the fastest
                        use of GCVSPL, since no iteration
                        is performed in p.
              |MD| = 2: Generalized cross validation.
              |MD| = 3: True predicted mean-squared error,
                        with prior given variance in VAL.
              |MD| = 4: Prior given number of degrees of
                        freedom in VAL (ZERO.le.VAL.le.N-M).
              After return from MD.ne.1, the same number of
              degrees of freedom can be obtained, for identical
              weight factors and knot positions, by selecting
              |MD|=1, and by copying the value of p from WK(4)
              into VAL. In this way, no iterative optimization
              is required when processing other data in Y.
      VAL     Mode value, as described above under MD.
      IDER    Derivative order required, with 0.le.IDER
              and IDER.le.2*M. If IDER.eq.0, the function
              value is returned; otherwise, the IDER-th
              derivative of the spline is returned.

Return values:
      YF(NF,1)Approximated values at XF.
      WK(IWK) On normal exit, the first 6 values of WK are
              assigned as follows:
              WK(1) = Generalized Cross Validation value
              WK(2) = Mean Squared Residual.
              WK(3) = Estimate of the number of degrees of
                      freedom of the residual sum of squares
                      per dataset, with 0.lt.WK(3).lt.N-M.
              WK(4) = Smoothing parameter p, multiplicative
                      with the splines' derivative constraint.
              WK(5) = Estimate of the true mean squared error
                      (different formula for |MD| = 3).
              WK(6) = Gauss-Markov error variance.

              If WK(4) -->  0 , WK(3) -->  0 , and an inter-
              polating spline is fitted to the data (p --> 0).
              A very small value > 0 is used for p, in order
              to avoid division by zero in the GCV function.

              If WK(4) --> inf, WK(3) --> N-M, and a least-
              squares polynomial of order M (degree M-1) is
              fitted to the data (p --> inf). For numerical
              reasons, a very high value is used for p.

              Upon return, the contents of WK can be used for
              covariance propagation in terms of the matrices
              B and WE: see the source listings. The variance
              estimate for dataset J follows as WK(6)/WY(J).

Remarks:
      (1) GCVSPL calculates a natural spline of order 2*M (degree
      2*M-1) which smoothes or interpolates a given set of data
      points, using statistical considerations to determine the
      amount of smoothing required (Craven & Wahba, 1979). If the
      error variance is a priori known, it should be supplied to
      the routine in VAL, for |MD|=3. The degree of smoothing is
      then determined to minimize an unbiased estimate of the true
      mean squared error. On the other hand, if the error variance
      is not known, one may select |MD|=2. The routine then deter-
      mines the degree of smoothing to minimize the generalized
      cross validation function. This is asymptotically the same
      as minimizing the true predicted mean squared error (Craven &
      Wahba, 1979). If the estimates from |MD|=2 or 3 do not appear
      suitable to the user (as apparent from the smoothness of the
      M-th derivative or from the effective number of degrees of
      freedom returned in WK(3) ), the user may select an other
      value for the noise variance if |MD|=3, or a reasonably large
      number of degrees of freedom if |MD|=4. If |MD|=1, the proce-
      dure is non-iterative, and returns a spline for the given
      value of the smoothing parameter p as entered in VAL.

      (2) The number of arithmetic operations and the amount of
      storage required are both proportional to N, so very large
      datasets may be accomodated. The data points do not have
      to be equidistant in the independant variable X or uniformly
      weighted in the dependant variable Y. However, the data
      points in X must be strictly increasing. Multiple dataset
      processing (K.gt.1) is numerically more efficient dan
      separate processing of the individual datasets (K.eq.1).

      (3) If |MD|=3 (a priori known noise variance), any value of
      N.ge.2*M is acceptable. However, it is advisable for N-2*M
      be rather large (at least 20) if |MD|=2 (GCV).

      (4) For |MD| > 1, GCVSPL tries to iteratively minimize the
      selected criterion function. This minimum is unique for |MD|
      = 4, but not necessarily for |MD| = 2 or 3. Consequently, 
      local optima rather that the global optimum might be found,
      and some actual findings suggest that local optima might
      yield more meaningful results than the global optimum if N
      is small. Therefore, the user has some control over the
      search procedure. If MD > 1, the iterative search starts
      from a value which yields a number of degrees of freedom
      which is approximately equal to N/2, until the first (local)
      minimum is found via a golden section search procedure
      (Utreras, 1980). If MD < -1, the value for p contained in
      WK(4) is used instead. Thus, if MD = 2 or 3 yield too noisy
      an estimate, the user might try |MD| = 1 or 4, for suitably
      selected values for p or for the number of degrees of
      freedom, and then run GCVSPL with MD = -2 or -3. The con-
      tents of N, M, K, X, WX, WY, and WK are assumed unchanged
      if MD < 0.

      (5) GCVSPL calculates the spline coefficient array C(N,K);
      this array can be used to calculate the spline function
      value and any of its derivatives up to the degree 2*M-1
      at any argument T within the knot range, using subrou-
      tines SPLDER and SEARCH, and the knot array X(N). Since
      the splines are constrained at their Mth derivative, only
      the lower spline derivatives will tend to be reliable
      estimates of the underlying, true signal derivatives.

      (6) GCVSPL combines elements of subroutine CRVO5 by Utre-
      ras (1980), subroutine SMOOTH by Lyche et al. (1983), and
      subroutine CUBGCV by Hutchinson (1985). The trace of the
      influence matrix is assessed in a similar way as described
      by Hutchinson & de Hoog (1985). The major difference is
      that the present approach utilizes non-symmetrical B-spline
      design matrices as described by Lyche et al. (1983); there-
      fore, the original algorithm by Erisman & Tinney (1975) has
      been used, rather than the symmetrical version adopted by
      Hutchinson & de Hoog.

References:
      P. Craven & G. Wahba (1979), Smoothing noisy data with
      spline functions. Numerische Mathematik 31, 377-403.

      A.M. Erisman & W.F. Tinney (1975), On computing certain
      elements of the inverse of a sparse matrix. Communications
      of the ACM 18(3), 177-179.

      M.F. Hutchinson & F.R. de Hoog (1985), Smoothing noisy data
      with spline functions. Numerische Mathematik 47(1), 99-106.

      M.F. Hutchinson (1985), Subroutine CUBGCV. CSIRO Division of
      Mathematics and Statistics, P.O. Box 1965, Canberra, ACT 2601,
      Australia.

      T. Lyche, L.L. Schumaker, & K. Sepehrnoori (1983), Fortran
      subroutines for computing smoothing and interpolating natural
      splines. Advances in Engineering Software 5(1), 2-5.

      F. Utreras (1980), Un paquete de programas para ajustar curvas
      mediante funciones spline. Informe Tecnico MA-80-B-209, Depar-
      tamento de Matematicas, Faculdad de Ciencias Fisicas y Matema-
      ticas, Universidad de Chile, Santiago.

      Wahba, G. (1980). Numerical and statistical methods for mildly,
      moderately and severely ill-posed problems with noisy data.
      Technical report nr. 595 (February 1980). Department of Statis-
      tics, University of Madison (WI), U.S.A.


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 35
B-spline data smoothing subroutine.





