Skip to content

Commit 3da8767

Browse files
JamesRamsay5JamesRamsay5
committed
update to fda_5.5.0
1 parent 54cc378 commit 3da8767

File tree

1 file changed

+268
-0
lines changed

1 file changed

+268
-0
lines changed

‎R/smooth.morph2.R

Lines changed: 268 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,268 @@
1+
smooth_morph2 <- function(x, y, ylim, WfdPar, wt=matrix(1,nobs,1),
2+
conv=1e-4, iterlim=20, dbglev=0) {
3+
#SMOOTH_MORPH smooths the relationship of Y to ARGVALS
4+
# by fitting a monotone fn. f(argvals) <- b_0 + b_1 D^{-1} exp W(t)
5+
# where W is a function defined over the same range as ARGVALS,
6+
# W + ln b_1 <- log Df and w <- D W <- D^2f/Df.
7+
# b_0 and b_1 are chosen so that values of f
8+
# are within the interval [ylim[1],ylim[2]].
9+
# The fitting criterion is penalized mean squared stop:
10+
# PENSSE(lambda) <- \sum [y_i - f(t_i)]^2 +
11+
# \lambda * \int [L W]^2
12+
# W(argvals) is expanded by the basis in functional data object Wfdobj.
13+
#
14+
# Arguments are
15+
# X argument value array
16+
# Y data array containing the the values to be fit
17+
# YLIM Ordinate value limits. The values of the estimated function
18+
# range between these limits. The abscissa range is
19+
# defined in WfdPar.
20+
# WFDPAR A functional parameter or fdPar object. This object
21+
# contains the specifications for the functional data
22+
# object to be estimated by smoothing the data. See
23+
# comment lines in function fdPar for details.
24+
# The functional data object WFD in FDPAROBJ is used
25+
# to initialize the optimization process.
26+
# It's coefficient array has a single column, and these
27+
# are the starting values for the iterative minimization
28+
# of mean squared stop.
29+
# This argument may also be either a FD object, or a
30+
# BASIS object. In this case, the smoothing parameter
31+
# LAMBDA is set to 0.
32+
# WT a vector of weights
33+
# CONV convergence criterion, 0.0001 by default
34+
# ITERLIM maximum number of iterations, 20 by default
35+
# DBGLEV Controls the level of output on each iteration. if (0,
36+
# no output, if (1, output at each iteration, if (higher,
37+
# output at each line search iteration. 1 by default.
38+
#
39+
# Returns are:
40+
# WFD Functional data object for W. It's coefficient vector
41+
# contains the optimized coefficients.
42+
43+
# last modified 26 October 2021
44+
45+
nobs <- length(x) # number of observations
46+
47+
# check consistency of x and y
48+
49+
if (length(y) != nobs) {
50+
stop('Arguments X and Y are not of same length')
51+
}
52+
x <- as.matrix(x)
53+
y <- as.matrix(y)
54+
55+
# check WfdPar
56+
57+
if (!is.fdPar(WfdPar)) {
58+
if (is.fd(WfdPar)) WfdPar <- fdPar(WfdPar)
59+
else stop(paste('WFDPAR is not a functional parameter object, ',
60+
'and not a functional data object.'))
61+
}
62+
63+
# extract information from WfdPar
64+
65+
Wfdobj <- WfdPar$fd
66+
Wbasis <- Wfdobj$basis # basis for Wfdobj
67+
Wnbasis <- Wbasis$nbasis # no. basis functions
68+
Wtype <- Wbasis$type
69+
xlim <- Wbasis$rangeval
70+
WLfdobj <- int2Lfd(WfdPar$Lfd)
71+
Wlambda <- WfdPar$lambda
72+
if (any(wt < 0)) {
73+
stop('One or more weights are negative.')
74+
}
75+
76+
cvec <- Wfdobj$coefs # initial coefficients
77+
78+
# first coefficient is not active for bspline or fourer bases
79+
80+
if (Wtype == 'bspline' || Wtype == 'fourier') {
81+
active <- 2:Wnbasis
82+
} else {
83+
active <- 1:Wnbasis
84+
}
85+
inactive <- rep(TRUE,Wnbasis)
86+
inactive[active] = FALSE
87+
88+
# check range of argvals
89+
90+
if (x[1] < xlim[1] || x[nobs] > xlim[2]) {
91+
stop('Values in ARGVALS are out of bounds.')
92+
}
93+
94+
# initialize matrix Kmat defining penalty term
95+
96+
if (Wlambda > 0) {
97+
Kmat <- Wlambda*eval.penalty(Wbasis, WLfdobj)
98+
} else {
99+
Kmat <- matrix(0,Wnbasis,Wnbasis)
100+
}
101+
102+
# load data into morphList
103+
104+
morphList <- list(x=x, y=y, xlim=xlim, ylim=ylim, wt=wt, inactive=inactive,
105+
Kmat=Kmat, Wlambda=Wlambda, Wbasis=Wbasis)
106+
107+
# Compute initial function and gradient values
108+
109+
result <- fngrad2(cvec, morphList)
110+
f <- result$PENSSE
111+
grad <- result$DPENSSE
112+
hmat <- result$D2PENSSE
113+
norm <- sqrt(sum(grad^2))
114+
115+
# compute initial badness of fit measures
116+
117+
fold <- f
118+
cvecold <- cvec
119+
120+
# evaluate the initial update vector
121+
122+
pvec <- -solve(hmat,grad)
123+
124+
# initialize iteration status arrays
125+
126+
iternum <- 0
127+
status <- c(iternum, fold, norm)
128+
if (dbglev >= 1) {
129+
cat("\nIter. PENSSE Grad Length")
130+
cat("\n")
131+
cat(iternum)
132+
cat(" ")
133+
cat(round(status[2],4))
134+
cat(" ")
135+
cat(round(status[3],4))
136+
}
137+
iterhist <- matrix(0,iterlim+1,length(status))
138+
iterhist[1,] <- status
139+
140+
# --------------------- Begin main iterations ---------------
141+
142+
STEPMAX <- 10
143+
itermax <- 20
144+
TOLX <- 1e-10
145+
fold <- f
146+
147+
# --------------- beginning of optimization loop -----------
148+
149+
for (iter in 1:iterlim) {
150+
iternum <- iternum + 1
151+
# line search
152+
result <- lnsrch(cvecold, fold, grad, pvec, fngrad2, morphList,
153+
STEPMAX, itermax, TOLX, dbglev)
154+
cvec <- result$x
155+
result <- fngrad2(cvec, morphList)
156+
f <- result$PENSSE
157+
grad <- result$DPENSSE
158+
hmat <- result$D2PENSSE
159+
norm <- sqrt(sum(grad^2))
160+
status <- c(iternum, f, norm)
161+
iterhist[iter+1,] <- status
162+
if (dbglev >= 1) {
163+
cat("\n")
164+
cat(iternum)
165+
cat(" ")
166+
cat(round(status[2],4))
167+
cat(" ")
168+
cat(round(status[3],4))
169+
}
170+
# test for convergence
171+
if (abs(f-fold) < conv) {
172+
break
173+
}
174+
if (f >= fold) {
175+
warning('Current function value does not decrease fit.')
176+
break
177+
}
178+
# evaluate the update vector
179+
pvec <- -solve(hmat,grad)
180+
cosangle <- -t(grad) %*% pvec/sqrt(sum(grad^2)*sum(pvec^2))
181+
if (cosangle < 0) {
182+
if (dbglev > 1) {
183+
print('cos(angle) negative')
184+
}
185+
pvec <- -grad
186+
}
187+
fold <- f
188+
cvecold <- cvec
189+
}
190+
191+
# construct output objects
192+
193+
Wfdobj <- fd(cvec, Wbasis)
194+
ywidth <- ylim[2] - ylim[1]
195+
hraw <- monfn( x, Wfdobj)
196+
hmax <- monfn( xlim[2], Wfdobj)
197+
hmin <- monfn( xlim[1], Wfdobj)
198+
hwidth <- hmax - hmin
199+
ysmth <- (xlim[1] - hmin) + ywidth*hraw/hwidth
200+
cat("\n")
201+
202+
return(list(Wfdobj=Wfdobj, f=f, grad=grad, hmat=hmat, norm=norm,
203+
ysmth=ysmth, iternum=iternum, iterhist=iterhist))
204+
}
205+
# ----------------------------------------------------------------
206+
207+
fngrad2 <- function(cvec, morphList) {
208+
# extract data from morphList
209+
x <- as.matrix(morphList$x)
210+
y <- as.matrix(morphList$y)
211+
xlim <- as.matrix(morphList$xlim)
212+
ylim <- as.matrix(morphList$ylim)
213+
wt <- as.matrix(morphList$wt)
214+
inactive <- as.matrix(morphList$inactive)
215+
Kmat <- morphList$Kmat
216+
Wlambda <- morphList$Wlambda
217+
Wbasis <- morphList$Wbasis
218+
219+
nobs <- length(x)
220+
ywidth <- ylim[2] - ylim[1]
221+
Wfdobj <- fd(cvec, Wbasis)
222+
Wnbasis <- Wbasis$nbasis
223+
# get unnormalized function and gradient values
224+
hraw <- as.matrix(monfn( x, Wfdobj))
225+
Dhraw <- mongrad(x, Wfdobj)
226+
# adjust functions and derivatives for normalization
227+
hmax <- monfn( xlim[2], Wfdobj)
228+
hmin <- monfn( xlim[1], Wfdobj)
229+
hwidth <- hmax - hmin
230+
h <- (xlim[1] - hmin) + ywidth*hraw/hwidth
231+
Dhmax <- mongrad(xlim[2], Wfdobj)
232+
Dhmin <- mongrad(xlim[1], Wfdobj)
233+
Dhwidth <- Dhmax - Dhmin
234+
Dh <- ywidth*(Dhraw*hwidth - hraw %*% Dhwidth)/hwidth^2
235+
res <- y - h
236+
f <- mean(res^2*wt)
237+
# if (required, gradient is computed
238+
temp <- Dh*(wt %*% matrix(1,1,Wnbasis))
239+
grad <- -2*t(temp) %*% res/nobs
240+
if (Wlambda > 0) {
241+
grad <- grad + 2 * Kmat %*% cvec
242+
f <- f + t(cvec) %*% Kmat %*% cvec
243+
}
244+
if (!is.null(inactive)) {
245+
grad[inactive] <- 0
246+
}
247+
# if (required, hessian matrix is computed
248+
wtroot <- sqrt(wt)
249+
temp <- Dh*(wtroot %*% matrix(1,1,Wnbasis))
250+
hmat <- 2*t(temp) %*% temp/nobs
251+
# adjust for penalty
252+
if (Wlambda > 0) {
253+
hmat <- hmat + 2*Kmat
254+
}
255+
# adjust for inactive coefficients
256+
if (!is.null(inactive)) {
257+
hmat[inactive, ] <- 0
258+
hmat[, inactive] <- 0
259+
hmat[inactive, inactive] <- diag(rep(1,Wnbasis)[inactive])
260+
}
261+
262+
PENSSE <- f
263+
DPENSSE <- grad
264+
D2PENSSE <- hmat
265+
266+
return(list(PENSSE=PENSSE, DPENSSE=DPENSSE, D2PENSSE=D2PENSSE))
267+
268+
}

0 commit comments

Comments
 (0)