-
Notifications
You must be signed in to change notification settings - Fork 5
/
drop.step_2d.R
76 lines (65 loc) · 2.57 KB
/
drop.step_2d.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
#' Function that tries dropping knots to find an improvement in fit
#'
#' @author Cameron Walker, Department of Engineering Science, University of Auckland.
#'
#' @export
#'
"drop.step_2d" <- function(radii,invInd,dists,explData,response,knotgrid,maxIterations,fitnessMeasure,
point,knotPoint,position,aR,BIC,track,out.lm,improveDrop,minKnots,tol=0,baseModel,radiusIndices,models, interactionTerm, data, initDisp, cv.opts, basis, printout) {
if (isS4(baseModel)){
attributes(baseModel@misc$formula)$.Environment<-environment()
} else {
attributes(baseModel$formula)$.Environment<-environment()
}
if(printout){
print("******************************************************************************")
print("Simplifying model...")
print("******************************************************************************")
}
fuse <- 0
improve <- 1
fitStat<-BIC[length(BIC)]
newRadii = radiusIndices
while ((improve) & (fuse < maxIterations) ) {
fuse <- fuse + 1
improve <- 0
if (length(aR) > minKnots) {
for (i in 1:length(aR)) {
tempR <- aR
tempR <- tempR[-i]
tempRadii = radiusIndices[-i]
output<-fit.thinPlate_2d(fitnessMeasure, dists,tempR,radii,baseModel,tempRadii,models, fitStat, interactionTerm, data, initDisp, cv.opts, basis)
initModel<-output$currentModel
models<-output$models
initBIC<-output$fitStat
out<-choose.radii(initBIC,1:length(tempRadii),tempRadii,radii,initModel,dists,tempR,baseModel,fitnessMeasure,response,models, interactionTerm, data, initDisp, cv.opts, basis, printout)
tempRadii=out$radiusIndices
tempOut.lm=out$out.lm
models=out$models
tempMeasure<-out$BIC
if (tempMeasure +tol < fitStat) {
out.lm <- tempOut.lm
fitStat<-tempMeasure
if(printout){
print("drop **********************************")
}
newR <- tempR
newRadii = tempRadii
tempKnot <- i
improve <- 1
improveDrop <- 1
}
}
if (improve) {
aR <- newR
point <- c(point,knotPoint[tempKnot])
position[knotPoint[tempKnot]]<-length(point)
knotPoint <- knotPoint[-tempKnot]
radiusIndices <- newRadii
}
}
}
### }
list(point=point,knotPoint=knotPoint,position=position,aR=aR,BIC=fitStat,track=track,out.lm=out.lm,improveDrop=improveDrop,
radiusIndices=radiusIndices,models=models)
}