Nothing
npscoefbw <-
function(...){
args = list(...)
if (is(args[[1]],"formula"))
UseMethod("npscoefbw",args[[1]])
else if (!is.null(args$formula))
UseMethod("npscoefbw",args$formula)
else
UseMethod("npscoefbw",args[[which(names(args)=="bws")[1]]])
}
npscoefbw.formula <-
function(formula, data, subset, na.action, call, ...){
orig.class <- if (missing(data))
sapply(eval(attr(terms(formula), "variables"), environment(formula)),class)
else sapply(eval(attr(terms(formula), "variables"), data, environment(formula)),class)
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "subset", "na.action"),
names(mf), nomatch = 0)
mf <- mf[c(1,m)]
if(!missing(call) && is.call(call)){
## rummage about in the call for the original formula
for(i in 1:length(call)){
if(tryCatch(class(eval(call[[i]])) == "formula",
error = function(e) FALSE))
break;
}
mf[[2]] <- call[[i]]
}
mf[[1]] <- as.name("model.frame")
chromoly <- explodePipe(mf[["formula"]])
bronze <- sapply(chromoly, paste, collapse = " + ")
mf[["formula"]] <-
as.formula(paste(bronze[1]," ~ ",
paste(bronze[2:length(bronze)],
collapse =" + ")),
env = environment(formula))
mf[["formula"]] <- terms(mf[["formula"]])
if(all(orig.class == "ts")){
args <- (as.list(attr(mf[["formula"]], "variables"))[-1])
attr(mf[["formula"]], "predvars") <- as.call(c(quote(as.data.frame),as.call(c(quote(ts.intersect), args))))
}else if(any(orig.class == "ts")){
arguments <- (as.list(attr(mf[["formula"]], "variables"))[-1])
arguments.normal <- arguments[which(orig.class != "ts")]
arguments.timeseries <- arguments[which(orig.class == "ts")]
ix <- sort(c(which(orig.class == "ts"),which(orig.class != "ts")),index.return = TRUE)$ix
attr(mf[["formula"]], "predvars") <- bquote(.(as.call(c(quote(cbind),as.call(c(quote(as.data.frame),as.call(c(quote(ts.intersect), arguments.timeseries)))),arguments.normal,check.rows = TRUE)))[,.(ix)])
}
mf <- eval(mf, parent.frame())
ydat <- model.response(mf)
xdat <- mf[, chromoly[[2]], drop = FALSE]
if (!(miss.z <- !(length(chromoly) == 3)))
zdat <- mf[, chromoly[[3]], drop = FALSE]
tbw <- eval(parse(text = paste('npscoefbw(xdat = xdat, ydat = ydat,',
ifelse(miss.z,'','zdat = zdat,'), '...)')))
## clean up (possible) inconsistencies due to recursion ...
tbw$call <- match.call(expand.dots = FALSE)
environment(tbw$call) <- parent.frame()
tbw$formula <- formula
tbw$rows.omit <- as.vector(attr(mf,"na.action"))
tbw$nobs.omit <- length(tbw$rows.omit)
tbw$terms <- attr(mf,"terms")
tbw$chromoly <- chromoly
tbw <-
updateBwNameMetadata(nameList =
list(ynames =
attr(mf, "names")[attr(tbw$terms, "response")]),
bws = tbw)
tbw
}
npscoefbw.NULL <-
function(xdat = stop("invoked without data 'xdat'"),
ydat = stop("invoked without data 'ydat'"),
zdat = NULL,
bws, ...){
miss.z <- missing(zdat)
xdat <- toFrame(xdat)
if(!miss.z)
zdat <- toFrame(zdat)
bws <- double(ifelse(miss.z, ncol(xdat), ncol(zdat)))
tbw <- eval(parse(text = paste("npscoefbw.default(xdat = xdat, ydat = ydat, bws = bws,",
ifelse(miss.z,"", "zdat = zdat,"), "...)")))
## clean up (possible) inconsistencies due to recursion ...
mc <- match.call(expand.dots = FALSE)
environment(mc) <- parent.frame()
tbw$call <- mc
tbw <-
updateBwNameMetadata(nameList = list(ynames = deparse(substitute(ydat))),
bws = tbw)
tbw
}
npscoefbw.scbandwidth <-
function(xdat = stop("invoked without data 'xdat'"),
ydat = stop("invoked without data 'ydat'"),
zdat = NULL,
bws,
nmulti,
random.seed = 42,
cv.iterate = FALSE,
cv.num.iterations = 1,
backfit.iterate = FALSE,
backfit.maxiter = 100,
backfit.tol = .Machine$double.eps,
bandwidth.compute = TRUE,
optim.method = c("Nelder-Mead", "BFGS", "CG"),
optim.maxattempts = 10,
optim.reltol = sqrt(.Machine$double.eps),
optim.abstol = .Machine$double.eps,
optim.maxit = 500,
...){
## Save seed prior to setting
if(exists(".Random.seed", .GlobalEnv)) {
save.seed <- get(".Random.seed", .GlobalEnv)
exists.seed = TRUE
} else {
exists.seed = FALSE
}
set.seed(random.seed)
miss.z <- missing(zdat)
xdat <- toFrame(xdat)
if (!miss.z)
zdat <- toFrame(zdat)
if (missing(nmulti)){
nmulti <- min(5,length(bws$bw))
}
if (!(is.vector(ydat) | is.factor(ydat)))
stop("'ydat' must be a vector or a factor")
eval(parse(text = paste("bwMatch(",
ifelse(miss.z, "xdat, bws$xdati", "zdat, bws$zdati"),")")))
if (dim(xdat)[1] != length(ydat))
stop("number of regression data and response data do not match")
if (ncol(xdat) == 1 && missing(cv.iterate))
cv.iterate = FALSE
if (cv.num.iterations < 1 & cv.iterate){
stop("invalid number of iterations specified")
}
if (!all(bws$xdati$icon))
stop("Only continuous 'x' regressors are supported in this version.")
optim.method <- match.arg(optim.method)
## catch and destroy NA's
goodrows = 1:dim(xdat)[1]
rows.omit =
eval(parse(text = paste("attr(na.omit(data.frame(xdat,ydat",
ifelse(miss.z,'',',zdat'),')), "na.action")')))
goodrows[rows.omit] = 0
if (all(goodrows==0))
stop("Data has no rows without NAs")
xdat = xdat[goodrows,,drop = FALSE]
ydat = ydat[goodrows]
if(!miss.z)
zdat <- zdat[goodrows,, drop = FALSE]
nrow = dim(xdat)[1]
ncol = dim(xdat)[2]
## at this stage, data to be sent to the c routines must be converted to
## numeric type.
if (is.factor(ydat))
ydat <- dlev(ydat)[as.integer(ydat)]
else
ydat <- as.double(ydat)
xdat <- toMatrix(xdat)
## if (!miss.z)
## zdat <- toMatrix(zdat)
## bad data
if (qr(xdat)$rank < ncol(xdat)){
stop("columns of the independent variable (xdat) are linearly dependent")
}
n <- nrow(xdat)
## ... do bandwidth selection
## construct 'W' matrix
## in the future one will be able to use a switch to npksum
## to emulate W
W <- as.matrix(data.frame(1,xdat))
yW <- as.matrix(data.frame(ydat,1,xdat))
if (miss.z){
zdat <- xdat
dati <- bws$xdati
}
else
dati <- bws$zdati
mysd <- EssDee(zdat[, dati$icon, drop = FALSE])
nconfac <- n^(-1.0/(2.0*bws$ckerorder+bws$ncon))
ncatfac <- n^(-2.0/(2.0*bws$ckerorder+bws$ncon))
bws$sdev <- mysd
bws$nconfac <- nconfac
bws$ncatfac <- ncatfac
total.time <-
system.time({
if (bandwidth.compute){
maxPenalty <- sqrt(.Machine$double.xmax)
overall.cv.ls <- function(param) {
sbw <- bws
sbw$bw <- param
sbw$bandwidth[[1]] <- param
if (!validateBandwidthTF(sbw) || ((bws$nord+bws$nuno > 0) && any(param[!bws$icon] > 2.0*x.scale[!bws$icon])))
return(maxPenalty)
bws$bw <- param
if(bws$scaling)
bws$bandwidth[[1]] <- sapply(1:bws$ndim, function(i) { bws$bw[i]*ifelse(bws$icon[i],nconfac*bws$sdev[sum(dati$icon[1:i])], ncatfac) })
else
bws$bandwidth[[1]] <- bws$bw
tww <- npksum(txdat = zdat, tydat = yW, weights = yW, bws = bws,
leave.one.out = TRUE)$ksum
mean.loo <- rep(maxPenalty,n)
epsilon <- 1.0/n
ridge <- double(n)
doridge <- !logical(n)
nc <- ncol(tww[-1,-1,1])
ridger <- function(i) {
doridge[i] <<- FALSE
ridge.val <- ridge[i]*tww[-1,1,i][1]/NZD(tww[-1,-1,i][1,1])
W[i,, drop = FALSE] %*% tryCatch(solve(tww[-1,-1,i]+diag(rep(ridge[i],nc)),
tww[-1,1,i]+c(ridge.val,rep(0,nc-1))),
error = function(e){
ridge[i] <<- ridge[i]+epsilon
doridge[i] <<- TRUE
return(rep(maxPenalty,nc))
})
}
while(any(doridge)){
iloo <- (1:n)[doridge]
mean.loo[iloo] <- sapply(iloo, ridger)
}
cv.console <<- printClear(cv.console)
##cv.console <<- printPush(msg = paste("param:", param), console = cv.console)
stopifnot(all(is.finite(mean.loo)))
if(!any(mean.loo == maxPenalty)){
fv <- sum((ydat-mean.loo)^2)/n
cv.console <<- printPush(msg = paste("fval:", signif(fv, digits = options('digits')$digits)), console = cv.console)
} else {
cv.console <<- printPush(msg = "near-singular system encountered, ridging", console = cv.console)
fv <- maxPenalty
}
return(ifelse(is.finite(fv),fv,maxPenalty))
}
scoef.looE <- parse(text = paste('npscoef(bws = bws, txdat = xdat, tydat = ydat,',
ifelse(miss.z,'','tzdat = zdat,'),
'leave.one.out = TRUE, iterate = TRUE,',
'maxiter = backfit.maxiter, tol = backfit.tol,',
'betas = TRUE)'))
partial.cv.ls <- function(param, partial.index) {
sbw <- bws
sbw$bw <- param
sbw$bandwidth[[1]] <- param
if (!validateBandwidthTF(sbw) || ((bws$nord+bws$nuno > 0) && any(param[!bws$icon] > 2.0*x.scale[!bws$icon])))
return(maxPenalty)
if (backfit.iterate){
bws$bw.fitted[,partial.index] <- param
scoef.loo <- eval(scoef.looE)
partial.loo <- W[,partial.index]*scoef.loo$beta[,partial.index]
} else {
bws$bw <- param
if(bws$scaling)
bws$bandwidth[[1]] <- sapply(1:bws$ndim, function(i) { bws$bw[i]*ifelse(bws$icon[i],nconfac*bws$sdev[sum(dati$icon[1:i])], ncatfac) })
else
bws$bandwidth[[1]] <- bws$bw
tww <- npksum(txdat=zdat,
tydat=cbind(partial.orig * W[,partial.index],W[,partial.index]^2),
weights=cbind(partial.orig * W[,partial.index],1),
bws=bws,
leave.one.out=TRUE)$ksum
partial.loo <- W[,partial.index]*tww[1,2,]/NZD(tww[2,2,])
}
fv <- sum((partial.orig - partial.loo)^2)/n
cv.console <<- printClear(cv.console)
cv.console <<- printPush(msg = paste("fval:",
signif(fv, digits = options('digits')$digits)),
console = cv.console)
return(ifelse(is.finite(fv),fv,maxPenalty))
}
## Now we implement multistarting
fval.min <- .Machine$double.xmax
numimp <- 0
value.overall <- numeric(nmulti)
x.scale <- sapply(1:bws$ndim, function(i){
if (dati$icon[i]){
return(1.059224*(ifelse(bws$scaling, 1.0, mysd[sum(dati$icon[1:i])]*nconfac)))
}
if (dati$iord[i])
return(0.5*oMaxL(dati$all.nlev[[i]], kertype = bws$okertype)*
ifelse(bws$scaling,ncatfac,1.0))
if (dati$iuno[i])
return(0.5*uMaxL(dati$all.nlev[[i]], kertype = bws$ukertype)*
ifelse(bws$scaling,ncatfac,1.0))
})
console <- newLineConsole()
optim.control <- list(abstol = optim.abstol,
reltol = optim.reltol,
maxit = optim.maxit)
for (i in 1:nmulti) {
console <- printPush(msg = paste(sep="", "Multistart ", i, " of ", nmulti, "... "), console)
cv.console <- newLineConsole(console)
if (i == 1) {
tbw <- x.scale
if(all(bws$bw != 0))
tbw <- bws$bw
} else {
tbw <- runif(bws$ndim, min=0.5, max=1.5)*x.scale
}
suppressWarnings(optim.return <- optim(tbw,
fn = overall.cv.ls,
control = optim.control))
attempts <- 0
while((optim.return$convergence != 0) && (attempts <= optim.maxattempts)) {
attempts <- attempts + 1
tbw <- runif(bws$ndim, min=0.5, max=1.5)*x.scale
optim.control <- lapply(optim.control, '*', 10.0)
suppressWarnings(optim.return <- optim(tbw,
fn = overall.cv.ls,
control = optim.control))
}
cv.console <- printClear(cv.console)
value.overall[i] <- optim.return$value
if(optim.return$value < fval.min) {
param <- optim.return$par
min.overall <- optim.return$value
fval.min <- min.overall ## Added by jracine Jul 22 2010
numimp.overall <- numimp + 1
best.overall <- i
}
if(i < nmulti)
console <- printPop(console)
else
console <- printClear(console)
}
param.overall <- bws$bw <- param
if(cv.iterate){
console <- newLineConsole()
n.part <- (ncol(xdat)+1)
bws$bw.fitted <- matrix(data = bws$bw, nrow = length(bws$bw), ncol = n.part)
## obtain matrix of alpha.hat | h0 and beta.hat | h0
scoef <- eval(parse(text = paste('npscoef(bws = bws, txdat = xdat, tydat = ydat,',
ifelse(miss.z, '', 'tzdat = zdat,'),
'iterate = FALSE, betas = TRUE)')))
resid.full <- ydat - scoef$mean
for(i in 1:cv.num.iterations){
console <- printPush(msg = paste(sep="", "backfitting iteration ", i, " of ", cv.num.iterations, "... "), console)
for(j in 1:n.part){
console <- printPush(msg = paste(sep="", "partial residual ", j, " of ", n.part, "... "), console)
cv.console <- newLineConsole(console)
## estimate partial residuals
partial.orig <- W[,j] * scoef$beta[,j] + resid.full
## minimise
suppressWarnings(optim.return <-
optim(tbw, fn = partial.cv.ls,
control = optim.control,
partial.index = j))
cv.console <- printClear(cv.console)
## grab parameter
bws$bw.fitted[,j] <- optim.return$par
if (backfit.iterate){
## re-estimate all betas
scoef <- eval(parse(text = paste('npscoef(bws = bws, txdat = xdat, tydat = ydat,',
ifelse(miss.z, '', 'tzdat = zdat,'),
'iterate = TRUE, maxiter = backfit.maxiter,',
'tol = backfit.tol, betas = TRUE)')))
resid.full <- ydat - scoef$mean
} else {
bws$bw <- bws$bw.fitted[,j]
## estimate new beta.hats
if(bws$scaling)
bws$bandwidth[[1]] <- sapply(1:bws$ndim, function(i) { bws$bw[i]*ifelse(bws$icon[i],nconfac*bws$sdev[sum(dati$icon[1:i])], ncatfac) })
else
bws$bandwidth[[1]] <- bws$bw
tww <- npksum(txdat=zdat,
tydat=cbind(partial.orig * W[,j],W[,j]^2),
weights=cbind(partial.orig * W[,j],1),
bws=bws)$ksum
scoef$beta[,j] <- tww[1,2,]/NZD(tww[2,2,])
bws$bw <- param.overall
## estimate new full residuals
resid.full <- partial.orig - W[,j] * scoef$beta[,j]
}
console <- printPop(console)
}
if(i < cv.num.iterations)
console <- printPop(console)
else
console <- printClear(console)
}
scoef.loo <- eval(parse(text = paste('npscoef(bws = bws, txdat = xdat, tydat = ydat,',
ifelse(miss.z,'', 'tzdat = zdat,'),
'iterate = TRUE, maxiter = backfit.maxiter,',
'tol = backfit.tol, leave.one.out = TRUE)$mean')))
bws$fval.fitted <- sum((ydat - scoef.loo)^2)/n
}
bws$fval = min.overall
bws$ifval = best.overall
bws$numimp = numimp.overall
bws$fval.vector = value.overall
}
})[1]
bws$sfactor <- bws$bandwidth <- bws$bw
nfactor <- nrow^(-2.0/(2.0*bws$ckerorder+bws$ncon))
dfactor <- EssDee(zdat[, dati$icon, drop = FALSE])*nrow^(-1.0/(2.0*bws$ckerorder+sum(dati$icon)))
if (bws$scaling) {
bws$bandwidth[dati$icon] <- bws$bandwidth[dati$icon]*dfactor
if(bws$nuno > 0)
bws$bandwidth[dati$iuno] <- bws$bandwidth[dati$iuno]*nfactor
if(bws$nord > 0)
bws$bandwidth[dati$iord] <- bws$bandwidth[dati$iord]*nfactor
} else {
bws$sfactor[dati$icon] <- bws$sfactor[dati$icon]/dfactor
if(bws$nuno > 0)
bws$sfactor[dati$iuno] <- bws$sfactor[dati$iuno]/nfactor
if(bws$nord > 0)
bws$sfactor[dati$iord] <- bws$sfactor[dati$iord]/nfactor
}
## Restore seed
if(exists.seed) assign(".Random.seed", save.seed, .GlobalEnv)
bws <- scbandwidth(bw = bws$bw,
bwmethod = bws$method,
bwscaling = bws$scaling,
bwtype = bws$type,
ckertype = bws$ckertype,
ckerorder = bws$ckerorder,
ukertype = bws$ukertype,
okertype = bws$okertype,
fval = bws$fval,
ifval = bws$ifval,
numimp = bws$numimp,
fval.vector = bws$fval.vector,
nobs = bws$nobs,
xdati = bws$xdati,
ydati = bws$ydati,
zdati = bws$zdati,
xnames = bws$xnames,
ynames = bws$ynames,
znames = bws$znames,
sfactor = bws$sfactor,
bandwidth = bws$bandwidth,
rows.omit = rows.omit,
bandwidth.compute = bandwidth.compute,
optim.method = optim.method,
total.time = total.time)
bws
}
npscoefbw.default <-
function(xdat = stop("invoked without data 'xdat'"),
ydat = stop("invoked without data 'ydat'"),
zdat = NULL,
bws,
nmulti,
random.seed,
cv.iterate,
cv.num.iterations,
backfit.iterate,
backfit.maxiter,
backfit.tol,
bandwidth.compute = TRUE,
## dummy arguments for scbandwidth()
bwmethod, bwscaling, bwtype,
ckertype, ckerorder,
ukertype, okertype,
optim.method, optim.maxattempts,
optim.reltol, optim.abstol, optim.maxit,
...){
miss.z <- missing(zdat)
xdat <- toFrame(xdat)
if(!(is.vector(ydat) | is.factor(ydat)))
stop("'ydat' must be a vector or a factor")
if(!miss.z)
zdat <- toFrame(zdat)
## first grab dummy args for scbandwidth() and perform 'bootstrap'
## bandwidth call
mc.names <- names(match.call(expand.dots = FALSE))
margs <- c("bwmethod", "bwscaling", "bwtype", "ckertype", "ckerorder",
"ukertype", "okertype")
m <- match(margs, mc.names, nomatch = 0)
any.m <- any(m != 0)
tbw <- eval(parse(text=paste("scbandwidth(bw = bws",
ifelse(any.m, ",",""),
paste(mc.names[m], ifelse(any.m,"=",""), mc.names[m], collapse=", "),
", nobs = dim(xdat)[1],",
"xdati = untangle(xdat),",
"ydati = untangle(data.frame(ydat)),",
"zdati = untangle(zdat),",
"xnames = names(xdat),",
"ynames = deparse(substitute(ydat)),",
"znames = names(zdat),",
"bandwidth.compute = bandwidth.compute)")))
## next grab dummies for actual bandwidth selection and perform call
mc.names <- names(match.call(expand.dots = FALSE))
margs <- c("zdat", "bandwidth.compute",
"nmulti",
"random.seed",
"cv.iterate",
"cv.num.iterations",
"backfit.iterate",
"backfit.maxiter",
"backfit.tol",
"optim.method", "optim.maxattempts",
"optim.reltol", "optim.abstol", "optim.maxit")
m <- match(margs, mc.names, nomatch = 0)
any.m <- any(m != 0)
tbw <- eval(parse(text=paste("npscoefbw.scbandwidth(xdat=xdat, ydat=ydat, bws=tbw",
ifelse(any.m, ",",""),
paste(mc.names[m], ifelse(any.m,"=",""), mc.names[m], collapse=", "),
")")))
mc <- match.call(expand.dots = FALSE)
environment(mc) <- parent.frame()
tbw$call <- mc
return(tbw)
}
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.