Nothing
npscoef <-
function(bws, ...){
args <- list(...)
if (!missing(bws)){
if (is.recursive(bws)){
if (!is.null(bws$formula) && is.null(args$txdat))
UseMethod("npscoef",bws$formula)
else if (!is.null(bws$call) && is.null(args$txdat))
UseMethod("npscoef",bws$call)
else if (!is.call(bws))
UseMethod("npscoef",bws)
else
UseMethod("npscoef",NULL)
} else {
UseMethod("npscoef", NULL)
}
} else {
UseMethod("npscoef", NULL)
}
}
npscoef.formula <-
function(bws, data = NULL, newdata = NULL, y.eval = FALSE, ...){
tt <- terms(bws)
m <- match(c("formula", "data", "subset", "na.action"),
names(bws$call), nomatch = 0)
tmf <- bws$call[c(1,m)]
tmf[[1]] <- as.name("model.frame")
tmf[["formula"]] <- tt
umf <- tmf <- eval(tmf, envir = environment(tt))
tydat <- model.response(tmf)
txdat <- tmf[, bws$chromoly[[2]], drop = FALSE]
if (!(miss.z <- !(length(bws$chromoly) == 3)))
tzdat <- tmf[, bws$chromoly[[3]], drop = FALSE]
if ((has.eval <- !is.null(newdata))) {
if (!y.eval){
tt <- delete.response(tt)
orig.class <- sapply(eval(attr(tt, "variables"), newdata, environment(tt)),class)
## delete.response clobbers predvars, which is used for timeseries objects
## so we need to reconstruct it
if(all(orig.class == "ts")){
args <- (as.list(attr(tt, "variables"))[-1])
attr(tt, "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(tt, "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(tt, "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)])
}else{
attr(tt, "predvars") <- attr(tt, "variables")
}
}
umf <- emf <- model.frame(tt, data = newdata)
if (y.eval)
eydat <- model.response(emf)
exdat <- emf[, bws$chromoly[[2]], drop = FALSE]
if (!miss.z)
ezdat <- emf[, bws$chromoly[[3]], drop = FALSE]
}
ev <-
eval(parse(text=paste("npscoef(txdat = txdat, tydat = tydat,",
ifelse(miss.z, '','tzdat = tzdat,'),
ifelse(has.eval,paste("exdat = exdat,",
ifelse(y.eval,"eydat = eydat,",""),
ifelse(miss.z,'', 'ezdat = ezdat,')),""),
"bws = bws, ...)")))
ev$omit <- attr(umf,"na.action")
ev$rows.omit <- as.vector(ev$omit)
ev$nobs.omit <- length(ev$rows.omit)
ev$mean <- napredict(ev$omit, ev$mean)
ev$merr <- napredict(ev$omit, ev$merr)
if(ev$residuals){
ev$resid <- naresid(ev$omit, ev$resid)
}
return(ev)
}
npscoef.call <-
function(bws, ...) {
eval(parse(text = paste('npscoef(txdat = eval(bws$call[["xdat"]], environment(bws$call)),',
'tydat = eval(bws$call[["ydat"]], environment(bws$call)),',
ifelse(!is.null(bws$zdati), 'tzdat = eval(bws$call[["zdat"]], environment(bws$call)),',''),
'bws = bws, ...)')))
}
npscoef.default <- function(bws, txdat, tydat, tzdat, ...) {
sc <- sys.call()
sc.names <- names(sc)
## here we check to see if the function was called with tdat =
## if it was, we need to catch that and map it to dat =
## otherwise the call is passed unadulterated to npudensbw
bws.named <- any(sc.names == "bws")
txdat.named <- any(sc.names == "txdat")
tydat.named <- any(sc.names == "tydat")
tzdat.named <- any(sc.names == "tzdat")
no.bws <- missing(bws)
no.txdat <- missing(txdat)
no.tydat <- missing(tydat)
no.tzdat <- missing(tzdat)
## if bws was passed in explicitly, do not compute bandwidths
if(txdat.named)
txdat <- toFrame(txdat)
## if(tydat.named)
## tydat <- toFrame(tydat)
if(tydat.named)
tzdat <- toFrame(tzdat)
sc.bw <- sc
sc.bw[[1]] <- quote(npscoefbw)
if(bws.named){
sc.bw$bandwidth.compute <- FALSE
}
ostxy <- c('txdat','tydat','tzdat')
nstxy <- c('xdat','ydat','zdat')
m.txy <- match(ostxy, names(sc.bw), nomatch = 0)
if(any(m.txy > 0)) {
names(sc.bw)[m.txy] <- nstxy[m.txy > 0]
}
tbw <- eval.parent(sc.bw)
## because of some ambiguities in how the function might be called
## we only drop up to two unnamed arguments, when sometimes dropping
## three would be appropriate. also, for simplicity, we don't allow
## for inconsistent mixes of named/unnamed arguments so bws is named
## or unnamed, and t[xyz]dat collectively either named or unnamed
tz.str <- ifelse(tzdat.named, ",tzdat = tzdat",
ifelse(no.tzdat, "", "tzdat"))
if(no.bws){
tx.str <- ",txdat = txdat"
ty.str <- ",tydat = tydat"
} else {
tx.str <- ifelse(txdat.named, ",txdat = txdat","")
ty.str <- ifelse(tydat.named, ",tydat = tydat","")
if((!bws.named) && (!txdat.named)){
ty.str <- ifelse(tydat.named, ",tydat = tydat",
ifelse(no.tydat, "", "tydat"))
}
}
eval(parse(text=paste("npscoef(bws = tbw", tx.str, ty.str, tz.str, ",...)")))
}
npscoef.scbandwidth <-
function(bws,
txdat = stop("training data 'txdat' missing"),
tydat = stop("training data 'tydat' missing"),
tzdat = NULL,
exdat,
eydat,
ezdat,
residuals = FALSE,
errors = TRUE,
iterate = TRUE,
maxiter = 100,
tol = .Machine$double.eps,
leave.one.out = FALSE,
betas = FALSE, ...){
miss.z <- missing(tzdat)
miss.ex = missing(exdat)
miss.ey = missing(eydat)
## if miss.ex then if !miss.ey then ey and tx must match, to get
## oos errors alternatively if miss.ey you get is errors if
## !miss.ex then if !miss.ey then ey and ex must match, to get oos
## errors alternatively if miss.ey you get NO errors since we
## don't evaluate on the training data
txdat <- toFrame(txdat)
if (!(is.vector(tydat) | is.factor(tydat)))
stop("'tydat' must be a vector or a factor")
if (!miss.z)
tzdat <- toFrame(tzdat)
if (!miss.ex){
exdat <- toFrame(exdat)
if (!miss.z)
ezdat <- toFrame(ezdat)
if (! txdat %~% exdat )
stop("'txdat' and 'exdat' are not similar data frames!")
if (!miss.ey){
if (dim(exdat)[1] != length(eydat))
stop("number of evaluation data 'exdat' and dependent data 'eydat' do not match")
}
} else if(!miss.ey) {
if (dim(txdat)[1] != length(eydat))
stop("number of training data 'txdat' and dependent data 'eydat' do not match")
}
if(iterate && !is.null(bws$bw.fitted) && !miss.ex){
warning("iteration is not supported for out of sample evaluations; using overall bandwidths")
iterate = FALSE
}
## catch and destroy NA's
goodrows = 1:dim(txdat)[1]
rows.omit =
eval(parse(text = paste("attr(na.omit(data.frame(txdat, tydat",
ifelse(miss.z,'',',tzdat'),')), "na.action")')))
goodrows[rows.omit] = 0
if (all(goodrows==0))
stop("Training data has no rows without NAs")
txdat = txdat[goodrows,,drop = FALSE]
tydat = tydat[goodrows]
if (!miss.z)
tzdat <- tzdat[goodrows,, drop = FALSE]
if (!miss.ex){
goodrows = 1:dim(exdat)[1]
rows.omit = eval(parse(text=paste('attr(na.omit(data.frame(exdat',
ifelse(miss.ey,"",",eydat"),
ifelse(miss.z, "",",ezdat"),
')), "na.action")')))
goodrows[rows.omit] = 0
exdat = exdat[goodrows,,drop = FALSE]
if (!miss.ey)
eydat = eydat[goodrows]
if (!miss.z)
ezdat <- ezdat[goodrows,, drop = FALSE]
if (all(goodrows==0))
stop("Evaluation data has no rows without NAs")
}
## convert tydat, eydat to numeric, from a factor with levels from the y-data
## used during bandwidth selection.
if (is.factor(tydat)){
tydat <- adjustLevels(as.data.frame(tydat), bws$ydati)[,1]
tydat <- (bws$ydati$all.dlev[[1]])[as.integer(tydat)]
}
else
tydat <- as.double(tydat)
if (miss.ey)
eydat <- double()
else {
if (is.factor(eydat)){
eydat <- adjustLevels(as.data.frame(eydat), bws$ydati)[,1]
eydat <- (bws$ydati$all.dlev[[1]])[as.integer(eydat)]
}
else
eydat <- as.double(eydat)
}
## re-assign levels in training and evaluation data to ensure correct
## conversion to numeric type.
txdat <- adjustLevels(txdat, bws$xdati)
if (!miss.z)
tzdat <- adjustLevels(tzdat, bws$zdati)
if (!miss.ex){
exdat <- adjustLevels(exdat, bws$xdati)
if (!miss.z)
ezdat <- adjustLevels(ezdat, bws$zdati)
}
## grab the evaluation data before it is converted to numeric
if(miss.ex){
teval <- txdat
if (!miss.z)
teval <- list(exdat = txdat, ezdat = tzdat)
} else {
teval <- exdat
if (!miss.z)
teval <- list(exdat = exdat, ezdat = ezdat)
}
## put the unordered, ordered, and continuous data in their own objects
## data that is not a factor is continuous.
txdat <- toMatrix(txdat)
if (!miss.ex){
exdat <- toMatrix(exdat)
}
## from this point on txdat and exdat have been recast as matrices
## construct 'W' matrix
W.train <- W <- as.matrix(data.frame(1,txdat))
yW <- as.matrix(data.frame(tydat,1,txdat))
if (miss.z){
tzdat <- txdat
if (!miss.ex)
ezdat <- exdat
}
tww <- eval(parse(text=paste("npksum(txdat = tzdat, tydat = yW, weights = yW,",
ifelse(miss.ex, "", "exdat = ezdat,"),
"bws = bws, leave.one.out = leave.one.out)$ksum")))
tyw <- tww[-1,1,,drop=FALSE]
dim(tyw) <- dim(tyw)[-2]
tww <- tww[-1,-1,,drop=FALSE]
tnrow <- nrow(txdat)
enrow <- ifelse(miss.ex, nrow(txdat), nrow(exdat))
if (!miss.ex)
W <- as.matrix(data.frame(1,exdat))
## ridging jracine Jan 28 2009
maxPenalty <- sqrt(.Machine$double.xmax)
coef.mat <- matrix(maxPenalty,ncol(W),enrow)
epsilon <- 1.0/enrow
ridge <- double(enrow)
doridge <- !logical(enrow)
nc <- ncol(tww[,,1])
ridger <- function(i) {
doridge[i] <<- FALSE
ridge.val <- ridge[i]*tyw[,i][1]/NZD(tww[,,i][1,1])
tryCatch(solve(tww[,,i]+diag(rep(ridge[i],nc)),
tyw[,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:enrow)[doridge]
coef.mat[,iloo] <- sapply(iloo, ridger)
}
if (do.iterate <- (iterate && !is.null(bws$bw.fitted) && miss.ex)){
resid <- tydat - sapply(1:enrow, function(i) { W[i,, drop = FALSE] %*% coef.mat[,i] })
i = 0
max.err <- .Machine$double.xmax
aydat <- abs(tydat) + .Machine$double.eps
n.part <- (ncol(txdat)+1)
while((max.err > tol) & ((i <- i + 1) <= maxiter)){
resid.old <- resid
for(j in 1:n.part){
## estimate partial residuals
partial <- W[,j] * coef.mat[j,] + resid
## use to calculate new beta implicitly
twww <- npksum(txdat=tzdat,
tydat=cbind(partial * W[,j],W[,j]^2),
weights=cbind(partial * W[,j],1),
bws=bws,
leave.one.out=leave.one.out)$ksum
coef.mat[j,] <- twww[1,2,]/NZD(twww[2,2,])
## estimate new full residuals
resid <- partial - W[,j] * coef.mat[j,]
## repeat for consistency ?
}
max.err <- max(abs(resid.old - resid)/aydat)
}
if (max.err > tol)
warning(paste("backfit iterations did not converge. max err= ", max.err,", tol= ", tol,", maxiter= ", maxiter, sep=''))
mean <- tydat - resid
} else {
mean <- sapply(1:enrow, function(i) { W[i,, drop = FALSE] %*% coef.mat[,i] })
}
if (!miss.ey) {
RSQ = RSQfunc(eydat, mean)
MSE = MSEfunc(eydat, mean)
MAE = MAEfunc(eydat, mean)
MAPE = MAPEfunc(eydat, mean)
CORR = CORRfunc(eydat, mean)
SIGN = SIGNfunc(eydat, mean)
} else if(miss.ex) {
RSQ = RSQfunc(tydat, mean)
MSE = MSEfunc(tydat, mean)
MAE = MAEfunc(tydat, mean)
MAPE = MAPEfunc(tydat, mean)
CORR = CORRfunc(tydat, mean)
SIGN = SIGNfunc(tydat, mean)
}
if(errors | (residuals & miss.ex)){
tywtm <- npksum(txdat = tzdat,
tydat = yW,
weights = yW,
bws = bws,
leave.one.out = leave.one.out)$ksum
tyw <- tywtm[-1,1,]
tm <- tywtm[-1,-1,]
mean.fit <- rep(maxPenalty,nrow(txdat))
epsilon <- 1.0/nrow(txdat)
ridge.tm <- double(nrow(txdat))
doridge <- !logical(nrow(txdat))
nc <- ncol(tm[,,1])
ridger <- function(i) {
doridge[i] <<- FALSE
ridge.val <- ridge.tm[i]*tyw[,i][1]/NZD(tm[,,i][1,1])
W.train[i,, drop = FALSE] %*% tryCatch(solve(tm[,,i]+diag(rep(ridge.tm[i],nc)),
tyw[,i]+c(ridge.val,rep(0,nc-1))),
error = function(e){
ridge.tm[i] <<- ridge.tm[i]+epsilon
doridge[i] <<- TRUE
return(rep(maxPenalty,nc))
})
}
while(any(doridge)){
ii <- (1:nrow(txdat))[doridge]
mean.fit[ii] <- sapply(ii, ridger)
}
u2.W <- (resid <- tydat - mean.fit)^2
}
if(errors){
## kernel^2 integrals
k <- (int.kernels[switch(bws$ckertype,
"truncated gaussian" = CKER_TGAUSS,
gaussian = CKER_GAUSS + bws$ckerorder/2 - 1,
epanechnikov = CKER_EPAN + bws$ckerorder/2 - 1,
uniform = CKER_UNI)+1])^length(bws$bw)
u2.W <- sapply(1:tnrow, function(i) { W.train[i,, drop=FALSE]*u2.W[i] })
u2.W <- t(u2.W)
V.hat <- eval(parse(text = paste("npksum(txdat = tzdat, tydat = W.train, weights = u2.W,",
ifelse(!miss.ex, "exdat = ezdat,", ""), "bws = bws, leave.one.out = leave.one.out)$ksum")))
## asymptotics rely on positive definite nature of tww (ie. M.eval) and V.hat
## so choleski decomposition is used to assure their veracity
merr <- sqrt(sapply(1:enrow, function(i){
cm <- chol2inv(chol(tww[,,i]+diag(rep(ridge[i],nc))))
k*(W[i,,drop=FALSE] %*% cm %*% V.hat[,,i] %*% cm %*% t(W[i,,drop=FALSE]))
}))
}
eval(parse(text=paste("smoothcoefficient(bws = bws, eval = teval",
", mean = mean,",
ifelse(errors & !do.iterate,"merr = merr,",""),
ifelse(betas, "beta = t(coef.mat),",""),
ifelse(residuals, "resid = resid,",""),
"residuals = residuals, betas = betas,",
"ntrain = nrow(txdat), trainiseval = miss.ex,",
ifelse(miss.ey && !miss.ex, "",
"xtra=c(RSQ,MSE,MAE,MAPE,CORR,SIGN)"),")")))
}
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.