1 |
x |
|
y |
|
est |
|
iter |
|
pyhat |
|
plotit |
|
fr |
|
xlab |
|
ylab |
|
zlab |
|
theta |
|
phi |
|
expand |
|
scale |
|
zscale |
|
xout |
|
eout |
|
outfun |
|
ticktype |
|
... |
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 77 78 | ##---- Should be DIRECTLY executable !! ----
##-- ==> Define data, use random,
##-- or do help(data=index) for the standard data sets.
## The function is currently defined as
function (x, y, est = tmean, iter = 10, pyhat = FALSE, plotit = TRUE,
fr = 1, xlab = "X", ylab = "Y", zlab = "", theta = 50, phi = 25,
expand = 0.5, scale = TRUE, zscale = TRUE, xout = FALSE,
eout = xout, outfun = out, ticktype = "simple", ...)
{
m <- elimna(cbind(x, y))
if (xout) {
flag <- outfun(x, plotit = FALSE)$keep
x = x[flag, ]
y = y[flag]
}
x <- as.matrix(x)
p <- ncol(x)
if (p == 1)
val <- rungen(x[, 1], y, est = est, pyhat = TRUE, plotit = plotit,
fr = fr, xlab = xlab, ylab = ylab, ...)$output
if (p > 1) {
library(MASS)
library(akima)
np <- p + 1
x <- m[, 1:p]
y <- m[, np]
fhat <- matrix(NA, ncol = p, nrow = length(y))
fhat.old <- matrix(NA, ncol = p, nrow = length(y))
res <- matrix(NA, ncol = np, nrow = length(y))
dif <- 1
for (i in 1:p) fhat.old[, i] <- rungenv2(x[, i], y, est = est,
pyhat = TRUE, plotit = FALSE, fr = fr, ...)
eval <- NA
for (it in 1:iter) {
for (ip in 1:p) {
res[, ip] <- y
for (ip2 in 1:p) {
if (ip2 != ip)
res[, ip] <- res[, ip] - fhat.old[, ip2]
}
fhat[, ip] <- rungenv2(x[, ip], res[, ip], est = est,
pyhat = TRUE, plotit = FALSE, fr = fr, ...)
}
eval[it] <- sum(abs(fhat/sqrt(sum(fhat^2)) - fhat.old/sqrt(sum(fhat.old^2))))
if (it > 1) {
itm <- it - 1
dif <- abs(eval[it] - eval[itm])
}
fhat.old <- fhat
if (dif < 0.01)
break
}
val <- apply(fhat, 1, sum)
aval <- est(y - val, ...)
val <- val + aval
if (plotit && p == 2) {
fitr <- val
iout <- c(1:length(fitr))
nm1 <- length(fitr) - 1
for (i in 1:nm1) {
ip1 <- i + 1
for (k in ip1:length(fitr)) if (sum(x[i, ] ==
x[k, ]) == 2)
iout[k] <- 0
}
fitr <- fitr[iout >= 1]
mkeep <- x[iout >= 1, ]
fitr <- interp(mkeep[, 1], mkeep[, 2], fitr)
persp(fitr, theta = theta, phi = phi, xlab = xlab,
ylab = ylab, zlab = zlab, expand = expand, scale = scale,
ticktype = ticktype)
}
}
if (!pyhat)
val <- "Done"
val
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.