Nothing
smooth.morph <- function(x, y, WfdPar, wt=rep(1,nobs),
conv=.0001, iterlim=20, dbglev=0)
{
#SMOOTH_MORPH smooths the relationship of Y to X
# by fitting a monotone fn. f(argvals) = b_0 + b_1 D^{-1} exp W(t)
# where W is a function defined over the same range as ARGVALS,
# W + ln b_1 = log Df and w = D W = D^2f/Df.
# b_0 and b_1 are chosen so that f(t_1) = y_1 and f(t_n) = y_n.
# The fitting criterion is penalized mean squared error:
# PENSSE(lambda) = \sum [y_i - f(t_i)]^2 +
# \lambda * \int [L W]^2
# W(x) is expanded by the basis in functional data object Wfd.
# Arguments:
# X ... vector of argument values
# Y ... vector of function values to be fit
# WFDPAR ... functional parameter object for W(x). The coefficient array
# for WFDPAROBJ$FD has a single column, and these are the starting
# values for the iterative minimization of mean squared error.
# WT ... a vector of weights
# CONV ... convergence criterion, 0.0001 by default
# ITERLIM ... maximum number of iterations, 20 by default
# DBGLEV ... Controls the level of output on each iteration. If 0,
# no output, if 1, output at each iteration, if higher, output
# at each line search iteration. 1 by default.
# Returns a list containing:
# WFDOBJ ... functional data object for W(x). Its coefficients are
# those that optimize fit.
# FLIST ... List containing final criterion value, gradient, and
# gradient norm.
# ITERNUM ... number of iterations
# ITERHIST ... ITERNUM+1 by 5 array containing iteration history
# last modified 13 may 2012 by Spencer Graves
# previously modified 2 May 2012 by Jim Ramsay
# initialize some arrays
x <- as.vector(x)
nobs <- length(x) # number of observations
# check WfdPar
if (!(inherits(WfdPar, "fdPar"))) stop(
"Argument WFDPAROBJ is not a functional data object.")
# extract information from WfdPar
Wfdobj <- WfdPar$fd
Lfdobj <- WfdPar$Lfd
lambda <- WfdPar$lambda
# check LFDOBJ
Lfdobj <- int2Lfd(Lfdobj)
# extract further information
basisobj <- Wfdobj$basis # basis for W(x)
nbasis <- basisobj$nbasis # number of basis functions
type <- basisobj$type
rangeval <- basisobj$rangeval
cvec <- Wfdobj$coefs
ncvec <- length(cvec)
if (type == "bspline" || type == "fourier") {
active <- c(FALSE, rep(TRUE, nbasis-1))
} else {
active <- rep(TRUE, nbasis)
}
# check some arguments
if (any(diff(x) <= 0)) stop("Arguments are not strictly increasing.")
if (x[1] < rangeval[1] | x[nobs] > rangeval[2]) stop(
"Argument values are out of range.")
if (any(wt < 0)) stop("One or more weights are negative.")
if (all(wt == 0)) stop("All weights are zero.")
# set up some variables
wtroot <- sqrt(wt)
climit <- c(-100*rep(1,nbasis), 100*rep(1,nbasis))
inact <- !active # indices of inactive coefficients
# set up cell for storing basis function values
JMAX <- 15
basislist <- vector("list", JMAX)
# initialize matrix Kmat defining penalty term
if (lambda > 0) {
Kmat <- lambda*getbasispenalty(basisobj, Lfdobj)
} else {
Kmat <- NULL
}
# Compute initial function and gradient values
fngradlist <- fngrad.smooth.morph(y, x, wt, Wfdobj, lambda,
Kmat, inact, basislist)
Flist <- fngradlist[[1]]
Dyhat <- fngradlist[[2]]
# compute the initial expected Hessian
hessmat <- hesscal.smooth.morph(Dyhat, wtroot, lambda, Kmat, inact)
# evaluate the initial update vector for correcting the initial cvec
result <- linesearch.smooth.morph(Flist, hessmat, dbglev)
deltac <- result[[1]]
cosangle <- result[[2]]
# initialize iteration status arrays
iternum <- 0
status <- c(iternum, Flist$f, Flist$norm)
if (dbglev >= 1) {
cat("\nIter. PENSSE Grad Length")
cat("\n")
cat(iternum)
cat(" ")
cat(round(status[2],4))
cat(" ")
cat(round(status[3],4))
}
# if (dbglev == 0 && iterlim > 1)
# cat("Progress: Each dot is an iteration\n")
iterhist <- matrix(0,iterlim+1,length(status))
iterhist[1,] <- status
if (iterlim == 0)
return ( list( "Wfdobj" = Wfdobj, "Flist" = Flist,
"iternum" = iternum, "iterhist" = iterhist ) )
# ------- Begin iterations -----------
MAXSTEPITER <- 10
MAXSTEP <- 100
trial <- 1
reset <- FALSE
linemat <- matrix(0,3,5)
cvecold <- cvec
Foldlist <- Flist
dbgwrd <- dbglev >= 2
for (iter in 1:iterlim)
{
# if (dbglev == 0 && iterlim > 1) cat(".")
iternum <- iternum + 1
# initialize logical variables controlling line search
dblwrd <- c(0,0)
limwrd <- c(0,0)
stpwrd <- 0
ind <- 0
ips <- 0
# compute slope at 0 for line search
linemat[2,1] <- sum(deltac*Flist$grad)
# normalize search direction vector
sdg <- sqrt(sum(deltac^2))
deltac <- deltac/sdg
linemat[2,1] <- linemat[2,1]/sdg
# initialize line search vectors
linemat[,1:4] <- outer(c(0, linemat[2,1], Flist$f),rep(1,4))
stepiter <- 0
if (dbglev >= 2) {
cat("\n")
cat(paste(" ", stepiter, " "))
cat(format(round(t(linemat[,1]),6)))
cat("\n")
}
# return with error condition if initial slope is nonnegative
if (linemat[2,1] >= 0) {
if (dbgwrd >= 2) print("Initial slope nonnegative.")
ind <- 3
break
}
# return successfully if initial slope is very small
if (linemat[2,1] >= -1e-7) {
if (dbglev >= 2) print("Initial slope too small")
ind <- 0
break
}
# first step set to trial
linemat[1,5] <- trial
# Main iteration loop for linesearch
cvecnew <- cvec
Wfdnew <- Wfdobj
for (stepiter in 1:MAXSTEPITER)
{
# ensure that step does not go beyond limits on parameters
limflg <- FALSE
# check the step size
result <-
stepchk(linemat[1,5], cvec, deltac, limwrd, ind,
climit, active, dbgwrd)
linemat[1,5] <- result[[1]]
ind <- result[[2]]
limwrd <- result[[3]]
if (linemat[1,5] <= 1e-7)
{
# Current step size too small ... terminate
if (dbglev >= 2) {
print("Stepsize too small")
# print(avec[5])
print(linemat[1, 5])
}
if (limflg) ind <- 1 else ind <- 4
break
}
# compute new function value and gradient
cvecnew <- cvec + linemat[1,5]*deltac
Wfdnew$coefs <- as.matrix(cvecnew)
fngradlist <- fngrad.smooth.morph(y, x, wt, Wfdnew, lambda,
Kmat, inact, basislist)
Flist <- fngradlist[[1]]
Dyhat <- fngradlist[[2]]
linemat[3,5] <- Flist$f
# compute new directional derivative
linemat[2,5] <- sum(deltac*Flist$grad)
if (dbglev >= 2) {
cat(paste(" ", stepiter, " "))
cat(format(round(t(linemat[,5]),6)))
cat("\n")
}
# compute next line search step, also test for convergence
result <- stepit(linemat, ips, dblwrd, MAXSTEP)
linemat <- result[[1]]
ips <- result[[2]]
ind <- result[[3]]
dblwrd <- result[[4]]
trial <- linemat[1,5]
# ind == 0 mean convergence
if (ind == 0 | ind == 5) break
}
# end iteration loop
cvec <- cvecnew
Wfdobj <- Wfdnew
# check that function value has not increased
if (Flist$f > Foldlist$f) {
# if it has, terminate iterations with a message
if (dbglev >= 2) {
cat("Criterion increased: ")
cat(format(round(c(Foldlist$f, Flist$f),4)))
cat("\n")
}
# reset parameters and fit
cvec <- cvecold
Wfdobj$coefs <- cvec
Flist <- Foldlist
deltac <- -Flist$grad
if (reset) {
# This is the second time in a row that
# this has happened ... quit
if (dbglev >= 2) cat("Reset twice, terminating.\n")
return ( list( "Wfdobj" = Wfdobj, "Flist" = Flist,
"iternum" = iternum, "iterhist" = iterhist ) )
} else {
reset <- TRUE
}
} else {
if (abs(Foldlist$f - Flist$f) < conv) {
# cat("\n")
break
}
cvecold <- cvec
Foldlist <- Flist
hessmat <- hesscal.smooth.morph(Dyhat, wtroot, lambda, Kmat, inact)
# udate the line search direction
result <- linesearch.smooth.morph(Flist, hessmat, dbglev)
deltac <- result[[1]]
cosangle <- result[[2]]
reset <- FALSE
}
# store iteration status
status <- c(iternum, Flist$f, Flist$norm)
iterhist[iter+1,] <- status
if (dbglev >= 1) {
cat("\n")
cat(iternum)
cat(" ")
cat(round(status[2],4))
cat(" ")
cat(round(status[3],4))
}
}
return ( list( "Wfdobj" = Wfdobj, "Flist" = Flist,
"iternum" = iternum, "iterhist" = iterhist ) )
}
# ----------------------------------------------------------------
linesearch.smooth.morph <- function(Flist, hessmat, dbglev)
{
deltac <- -symsolve(hessmat,Flist$grad)
cosangle <- -sum(Flist$grad*deltac)/sqrt(sum(Flist$grad^2)*sum(deltac^2))
if (dbglev >= 2) {
cat(paste("\nCos(angle) =",format(round(cosangle,4))))
if (cosangle < 1e-7) {
if (dbglev >=2) cat("\nCosine of angle too small\n")
deltac <- -Flist$grad
}
}
return(list(deltac, cosangle))
}
# ----------------------------------------------------------------
fngrad.smooth.morph <- function(y, x, wt, Wfdobj, lambda,
Kmat, inact, basislist)
{
nobs <- length(x)
width <- y[nobs] - y[1]
cvec <- Wfdobj$coefs
nbasis <- length(cvec)
h <- monfn(x, Wfdobj)
Dyhat <- mongrad(x, Wfdobj)
# adjust h and Dyhat for normalization
hmax <- h[nobs]
Dymax <- Dyhat[nobs,]
Dyhat <- width*(hmax*Dyhat - outer(h, Dymax))/hmax^2
h <- y[1] + width*h/hmax
# update residuals and function values
res <- y - h
f <- mean(res^2*wt)
temp <- Dyhat*outer(wt,rep(1,nbasis))
grad <- -2*crossprod(temp,res)/nobs
if (lambda > 0) {
grad <- grad + 2 * Kmat %*% cvec
f <- f + t(cvec) %*% Kmat %*% cvec
}
if (any(inact)) grad[inact] <- 0
norm <- sqrt(sum(grad^2)) # gradient norm
Flist <- list("f"=f,"grad"=grad,"norm"=norm)
return(list(Flist, Dyhat))
}
# ----------------------------------------------------------------
hesscal.smooth.morph <- function(Dyhat, wtroot, lambda, Kmat, inact)
{
Dydim <- dim(Dyhat)
nobs <- Dydim[1]
nbasis <- Dydim[2]
temp <- Dyhat*(wtroot %*% matrix(1,1,nbasis))
hessmat <- 2*crossprod(temp)/nobs
# adjust for penalty
if (lambda > 0) hessmat <- hessmat + 2*Kmat
# adjust for inactive coefficients
if (any(inact)) {
eyemat <- diag(rep(1,nbasis))
hessmat[inact, ] <- 0
hessmat[ ,inact] <- 0
hessmat[inact,inact] <- eyemat[inact,inact]
}
return(hessmat)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.