dcv Function

public function dcv(xval, ndata, nconst, nord, nbkpt, bkpt, w)

dcv is a companion function subprogram for dfc. The documentation for dfc has complete usage instructions.

dcv is used to evaluate the variance function of the curve obtained by the constrained B-spline fitting subprogram, dfc. The variance function defines the square of the probable error of the fitted curve at any point, XVAL. One can use the square root of this variance function to determine a probable error band around the fitted curve.

dcv is used after a call to dfc. MODE, an input variable to dfc, is used to indicate if the variance function is desired. In order to use dcv, MODE must equal 2 or 4 on input to dfc. MODE is also used as an output flag from dfc. Check to make sure that MODE = 0 after calling dfc, indicating a successful constrained curve fit. The array SDDATA, as input to dfc, must also be defined with the standard deviation or uncertainty of the Y values to use dcv.

To evaluate the variance function after calling dfc as stated above, use dcv as shown here

VAR = DCV(XVAL,NDATA,NCONST,NORD,NBKPT,BKPT,W)

The variance function is given by

VAR = (transpose of B(XVAL))*C*B(XVAL)/DBLE(MAX(NDATA-N,1))

where N = NBKPT - NORD.

The vector B(XVAL) is the B-spline basis function values at X=XVAL. The covariance matrix, C, of the solution coefficients accounts only for the least squares equations and the explicitly stated equality constraints. This fact must be considered when interpreting the variance function from a data fitting problem that has inequality constraints on the fitted curve.

All the variables in the calling sequence for dcv are used in dfc except the variable XVAL. Do not change the values of these variables between the call to dfc and the use of dcv.

Reference

  • R. J. Hanson, Constrained least squares curve fitting to discrete data using B-splines, a users guide, Report SAND78-1291, Sandia Laboratories, December 1978.

Revision history

  • 780801 DATE WRITTEN. Hanson, R. J., (SNLA)
  • 890531 Changed all specific intrinsics to generic. (WRB)
  • 890831 Modified array declarations. (WRB)
  • 890911 Removed unnecessary intrinsics. (WRB)
  • 891006 Cosmetic changes to prologue. (WRB)
  • 891006 REVISION DATE from Version 3.2
  • 891214 Prologue converted to Version 4.0 format. (BAB)
  • 920501 Reformatted the REFERENCES section. (WRB)

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: xval

The point where the variance is desired
integer, intent(in) :: ndata

The number of discrete (X,Y) pairs for which dfc calculated a piece-wise polynomial curve.
integer, intent(in) :: nconst

The number of conditions that constrained the B-spline in dfc.
integer, intent(in) :: nord

The order of the B-spline used in dfc. The value of NORD must satisfy 1 < NORD < 20 .

(The order of the spline is one more than the degree of the piece-wise polynomial defined on each interval. This is consistent with the B-spline package convention. For example, NORD=4 when we are using piece-wise cubics.)
integer, intent(in) :: nbkpt

The number of knots in the array BKPT(). The value of NBKPT must satisfy NBKPT .GE. 2NORD.
real(kind=wp), intent(in) :: bkpt(*)

The array of knots. Normally the problem data interval will be included between the limits BKPT(NORD) and BKPT(NBKPT-NORD+1). The additional end knots BKPT(I),I=1,...,NORD-1 and I=NBKPT-NORD+2,...,NBKPT, are required by dfc to compute the functions used to fit the data.
real(kind=wp) :: w(*)

Real work array as used in dfc. See dfc for the required length of W(). The contents of W() must not be modified by the user if the variance function is desired.

Return Value real(kind=wp)


Calls

proc~~dcv~~CallsGraph proc~dcv bspline_defc_module::dcv proc~ddot bspline_blas_module::ddot proc~dcv->proc~ddot proc~dfspvn bspline_defc_module::dfspvn proc~dcv->proc~dfspvn

Source Code

   real(wp) function dcv (xval, ndata, nconst, nord, nbkpt, bkpt, w)

      real(wp),intent(in) :: xval !! The point where the variance is desired
      integer,intent(in) :: nbkpt !! The number of knots in the array BKPT(*).
                                  !! The value of NBKPT must satisfy NBKPT .GE. 2*NORD.
      integer,intent(in) :: nconst !! The number of conditions that constrained the B-spline in
                                   !! [[dfc]].
      integer,intent(in) :: ndata !! The number of discrete (X,Y) pairs for which [[dfc]]
                                  !! calculated a piece-wise polynomial curve.
      integer,intent(in) :: nord !! The order of the B-spline used in [[dfc]].
                                 !! The value of NORD must satisfy 1 < NORD < 20 .
                                 !!
                                 !! (The order of the spline is one more than the degree of
                                 !! the piece-wise polynomial defined on each interval.  This
                                 !! is consistent with the B-spline package convention.  For
                                 !! example, NORD=4 when we are using piece-wise cubics.)
      real(wp),intent(in) :: bkpt(*) !! The array of knots.  Normally the problem
                                     !! data interval will be included between the limits
                                     !! BKPT(NORD) and BKPT(NBKPT-NORD+1).  The additional end
                                     !! knots BKPT(I),I=1,...,NORD-1 and I=NBKPT-NORD+2,...,NBKPT,
                                     !! are required by [[dfc]] to compute the functions used to
                                     !! fit the data.
      real(wp) :: w(*) !! Real work array as used in [[dfc]].  See [[dfc]]
                       !! for the required length of W(*).  The contents of W(*)
                       !! must not be modified by the user if the variance function
                       !! is desired.

      real(wp) :: v(40)
      integer :: i, ileft, ip, is, last, mdg, mdw, n
      integer :: dfspvn_j
      real(wp), dimension(20) :: dfspvn_deltam, dfspvn_deltap

      real(wp),parameter :: zero = 0.0_wp

      ! set up variables for dfspvn
      dfspvn_j = 1
      dfspvn_deltam = zero
      dfspvn_deltap = zero

      mdg = nbkpt - nord + 3
      mdw = nbkpt - nord + 1 + nconst
      is = mdg*(nord + 1) + 2*max(ndata,nbkpt) + nbkpt + nord**2
      last = nbkpt - nord + 1
      ileft = nord
      do
         if (xval < bkpt(ileft+1) .or. ileft >= last - 1) exit
         ileft = ileft + 1
      end do
      call dfspvn(bkpt,nord,1,xval,ileft,v(nord+1), &
                  dfspvn_j, dfspvn_deltam, dfspvn_deltap)
      ileft = ileft - nord + 1
      ip = mdw*(ileft - 1) + ileft + is
      n = nbkpt - nord
      do i = 1, nord
         v(i) = ddot(nord,w(ip),1,v(nord+1),1)
         ip = ip + mdw
      end do
      dcv = max(ddot(nord,v,1,v(nord+1),1),zero)

      ! scale the variance so it is an unbiased estimate.
      dcv = dcv/max(ndata-n,1)

   end function dcv