R/df2lambda.R

In fda: Functional Data Analysis

df2lambda <- function(argvals, basisobj, wtvec=rep(1,n), Lfdobj=0, df=nbasis)
{
#  Convert a degree of freedom DF for a smooth to the equivalent value
#    of the smoothing parameter lambda.

#  Last modified 26 October 2005

n <- length(argvals)
nbasis <- basisobj$nbasis
if (df >= nbasis) {
   lambda <- 0
   return(lambda)
}

TOL    <- 1e-3
GOLD   <- 1.0
GLIMIT <- 2.0
TINY   <- 1.0e-20

#  find machine precision
eps <- 1
tol1 <- 1 + eps
while (tol1 > 1) {
   eps  <- eps/2
   tol1 <- 1 + eps
}
eps <- sqrt(eps)

#  ------  initialization of lambda by finding bracketing values ------------
#             a < b < c such that  fb < fa  and  fb < fc
#  first use input value for lambda unless it is zero, in which case -1
bx <- -4.0
#  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
lambda <- 10^(bx)
fb <- (lambda2df(argvals, basisobj, wtvec, Lfdobj, lambda) - df)^2
#  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#  now try bracketing the minimum by using a large value and a small
#  value.  If (this doesn't work, revert to the iterative method
#  at statement 5
if (bx >= -10 &&  bx <= 5) {
#  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   cx <- 5  #  the upper limit
#  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   lambda <- 10^(cx)
   fc <- (lambda2df(argvals, basisobj, wtvec, Lfdobj, lambda) - df)^2
#  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   ax <- -8  #  the lower limit
#  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   lambda <- 10^(ax)
   fa <- (lambda2df(argvals, basisobj, wtvec, Lfdobj, lambda) - df)^2
}
#  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#  check to see if minimum bracketed
#print(c(ax,bx,cx,fa,fb,fc, lambda))
if (fb >= fa || fb >= fc) {
  #  Failure to bracket minimum, proceed with iterative search for
  #    bracketing values.
  #  First, as an alternative value for ax, use the input value plus 0.1
  ax <- bx + 1
  #  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  lambda <- 10^(ax)
  fa <- (lambda2df(argvals, basisobj, wtvec, Lfdobj, lambda) - df)^2
  #  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #  now the bracketing process begins
  if (fb > fa) {
     #  exchange ax and bx
     dum <- ax
     ax  <- bx
     bx  <- dum
     dum <- fb
     fb  <- fa
     fa  <- dum
  }
  #  first guess at cx
  cx <- bx + GOLD*(bx - ax)
  #  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  lambda <- 10^(cx)
  fc <- (lambda2df(argvals, basisobj, wtvec, Lfdobj, lambda) - df)^2
  #  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #  check if (three values bracket minimum
  #print(c(ax,bx,cx,fa,fb,fc, lambda))
  while (fb >= fc) {
     r <- (bx - ax)*(fb - fc)
     q <- (bx - cx)*(fb - fa)
     u <- bx -
       ((bx - cx)*q - (bx - ax)*r)/(2.0*sign(max(c(abs(q-r),TINY)))*(q-r))
     ulim <- bx + GLIMIT*(cx - bx)
     if ((bx-u)*(u-cx) > 0.0) {
        #  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        lambda <- 10^(u)
        fu <- (lambda2df(argvals, basisobj, wtvec, Lfdobj, lambda) - df)^2
        #  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        if (fu < fc) {
           #  success
           ax <- bx
           bx <- u
           fa <- fb
           fb <- fu
           break
        }
        if (fu > fb) {
           #  also success
           cx <- u
           fc <- fu
           break
        }
        #  failure:  fu >= fb
        u <- cx + GOLD*(cx - bx)
        #  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        lambda <- 10^(u)
        fu <- (lambda2df(argvals, basisobj, wtvec, Lfdobj, lambda) - df)^2
        #  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     }
     if ((cx - u)*(u - ulim) > 0.0) {
        #  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        lambda <- 10^(u)
        fu <- (lambda2df(argvals, basisobj, wtvec, Lfdobj, lambda) - df)^2
        #  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        if (fu < fc) {
           bx <- cx
           cx <-  u
           u  <- cx + GOLD*(cx - bx)
           fb <- fc
           fc <- fu
           #  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
           lambda <- 10^(u)
           fu <- (lambda2df(argvals, basisobj, wtvec, Lfdobj, lambda) - df)^2
           #  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        }
     }
     if ((u-ulim)*(ulim-cx) >= 0.0) {
        u <- ulim
        #  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        lambda <- 10^(u)
        fu <- (lambda2df(argvals, basisobj, wtvec, Lfdobj, lambda) - df)^2
        #  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     } else {
        u <- cx + GOLD*(cx - bx)
        #  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        lambda <- 10^(u)
        fu <- (lambda2df(argvals, basisobj, wtvec, Lfdobj, lambda) - df)^2
        #  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     }
     ax <- bx
     bx <- cx
     cx <- u
     fa <- fb
     fb <- fc
     fc <- fu
     #print(c(ax,bx,cx,fa,fb,fc, lambda))
  }  #  end of while loop
}
#  ---------------------------------------------------------------------
#  --------------------  bracketing successful  ------------------------
#  ---------------------------------------------------------------------
a  <- min(c(ax,cx))
b  <- max(c(ax,cx))
v  <- bx
w  <- v
x  <- v
e  <- 0.0
fx <- fb
fv <- fx
fw <- fx
#  ---------------------------------------------------------------------
#  --------------------  main loop starts here -------------------------
#  ---------------------------------------------------------------------
xm   <- 0.5*(a + b)
tol1 <- eps*abs(x) + TOL/3
tol2 <- 2*tol1
crit <- abs(x - xm) - (tol2 - 0.5*(b - a))
#print(c(crit, lambda))
while (crit > 0) {
   #  is golden-section necessary?
   if (abs(e) > tol1) {
      #  fit parabola
      r <- (x - w)*(fx - fv)
      q <- (x - v)*(fx - fw)
      p <- (x - v)*q - (x - w)*r
      q <- 2.0*(q - r)
      if (q > 0.0)   p <- -p
      q <- abs(q)
      s <- e
      e <- d
      #  is parabola acceptable?
      if (abs(p) < abs(0.5*q*s) & p > q*(a - x) & p < q*(b - x)) {
         #  a parabolic interpolation step
         d <- p/q
         u <- x + d
         #  f must not be evaluated too close to a or b
         if ((u - a) < tol2 ||  b - u < tol2) {
            if (xm - x >= 0.0) d <- tol1 else d <- -tol1
         }
      } else {
         #  a golden-section step
         if (x >= xm) e <- a - x else e <- b - x
         d <- 0.382*e
      }
   } else {
      #  a golden-section step
      if (x >= xm) e <- a - x else e <- b - x
      d <- 0.382*e
   }
#  f must not be evaluated too close to x
   if (abs(d) >=  tol1) {
      u <- x + d
   } else {
      if (d >= 0.0) u <- x + tol1 else u <- x - tol1
  }
  #  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  lambda <- 10^u
  fu <- (lambda2df(argvals, basisobj, wtvec, Lfdobj, lambda) - df)^2
  #  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #  update  a, b, v, w, and x
  if (fu <= fx) {
     if (u  >= x) a <- x else b <- x
     v  <- w
     w  <- x
     x  <- u
     fv <- fw
     fw <- fx
     fx <- fu
  } else {
     if (u  < x) a <- u else b <- u
     if (fu <= fw || w == x) {
        v  <- w
        w  <- u
        fv <- fw
        fw <- fu
     }
     if (fu <= fv || v == x || v == w) {
        v  <- u
        fv <- fu
     }
  }
  xm   <- 0.5*(a + b)
  tol1 <- eps*abs(x) + TOL/3
  tol2 <- 2*tol1
  crit <- abs(x - xm) - (tol2 - 0.5*(b - a))
  #print(c(crit, lambda))
#  -------------------  } of main loop  ------------------------------
}
return(lambda)
}