[Returnanalytics-commits] r2713 - in pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7: . Code Literature
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Aug 3 12:10:41 CEST 2013
Author: shubhanm
Date: 2013-08-03 12:10:40 +0200 (Sat, 03 Aug 2013)
New Revision: 2713
Added:
pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/glmi.R
pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/nlsi.R
pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Literature/Thumbs.db
pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/lmi.R
Log:
Week 6-7 : Regression methods support for HC and HAc for lm(),glm, and nls.
Status & Pending work : Still some errors for nls and vignette to be added.Also comparitive test with matlab equivalent function
Added: pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/glmi.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/glmi.R (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/glmi.R 2013-08-03 10:10:40 UTC (rev 2713)
@@ -0,0 +1,92 @@
+glmi <- function (formula, family = gaussian, data,vcov = NULL, weights, subset,
+ na.action, start = NULL, etastart, mustart, offset, control = list(...),
+ model = TRUE, method = "glm.fit", x = FALSE, y = TRUE, contrasts = NULL,
+ ...)
+{
+ call <- match.call()
+ if (is.character(family))
+ family <- get(family, mode = "function", envir = parent.frame())
+ if (is.function(family))
+ family <- family()
+ if (is.null(family$family)) {
+ print(family)
+ stop("'family' not recognized")
+ }
+ if (missing(data))
+ data <- environment(formula)
+ mf <- match.call(expand.dots = FALSE)
+ m <- match(c("formula", "data", "subset", "weights", "na.action",
+ "etastart", "mustart", "offset"), names(mf), 0L)
+ mf <- mf[c(1L, m)]
+ mf$drop.unused.levels <- TRUE
+ mf[[1L]] <- as.name("model.frame")
+ mf <- eval(mf, parent.frame())
+ if (identical(method, "model.frame"))
+ return(mf)
+ if (!is.character(method) && !is.function(method))
+ stop("invalid 'method' argument")
+ if (identical(method, "glm.fit"))
+ control <- do.call("glm.control", control)
+ mt <- attr(mf, "terms")
+ Y <- model.response(mf, "any")
+ if (length(dim(Y)) == 1L) {
+ nm <- rownames(Y)
+ dim(Y) <- NULL
+ if (!is.null(nm))
+ names(Y) <- nm
+ }
+ X <- if (!is.empty.model(mt))
+ model.matrix(mt, mf, contrasts)
+ else matrix(, NROW(Y), 0L)
+ weights <- as.vector(model.weights(mf))
+ if (!is.null(weights) && !is.numeric(weights))
+ stop("'weights' must be a numeric vector")
+ if (!is.null(weights) && any(weights < 0))
+ stop("negative weights not allowed")
+ offset <- as.vector(model.offset(mf))
+ if (!is.null(offset)) {
+ if (length(offset) != NROW(Y))
+ stop(gettextf("number of offsets is %d should equal %d (number of observations)",
+ length(offset), NROW(Y)), domain = NA)
+ }
+ mustart <- model.extract(mf, "mustart")
+ etastart <- model.extract(mf, "etastart")
+ fit <- eval(call(if (is.function(method)) "method" else method,
+ x = X, y = Y, weights = weights, start = start, etastart = etastart,
+ mustart = mustart, offset = offset, family = family,
+ control = control, intercept = attr(mt, "intercept") >
+ 0L))
+ if (length(offset) && attr(mt, "intercept") > 0L) {
+ fit2 <- eval(call(if (is.function(method)) "method" else method,
+ x = X[, "(Intercept)", drop = FALSE], y = Y, weights = weights,
+ offset = offset, family = family, control = control,
+ intercept = TRUE))
+ if (!fit2$converged)
+ warning("fitting to calculate the null deviance did not converge -- increase 'maxit'?")
+ fit$null.deviance <- fit2$deviance
+ }
+ if (model)
+ fit$model <- mf
+ fit$na.action <- attr(mf, "na.action")
+ if (x)
+ fit$x <- X
+ if (!y)
+ fit$y <- NULL
+ fit <- c(fit, list(call = call, formula = formula, terms = mt,
+ data = data, offset = offset, control = control, method = method,
+ contrasts = attr(X, "contrasts"), xlevels = .getXlevels(mt,
+ mf)))
+ class(fit) <- c(fit$class, c("glm", "lm"))
+ fit
+ if(is.null(vcov)) {
+ se <- vcov(fit)
+ } else {
+ if (is.function(vcov))
+ se <- vcov(fit)
+ else
+ se <- vcov
+ }
+ fit = list(fit,vHaC = se)
+ fit
+
+}
\ No newline at end of file
Added: pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/nlsi.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/nlsi.R (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/nlsi.R 2013-08-03 10:10:40 UTC (rev 2713)
@@ -0,0 +1,178 @@
+nlsi <- function (formula, data = parent.frame(),vcov = NULL, start, control = nls.control(),
+ algorithm = c("default", "plinear", "port"), trace = FALSE,
+ subset, weights, na.action, model = FALSE, lower = -Inf,
+ upper = Inf, ...)
+{
+ formula <- as.formula(formula)
+ algorithm <- match.arg(algorithm)
+ if (!is.list(data) && !is.environment(data))
+ stop("'data' must be a list or an environment")
+ mf <- match.call()
+ varNames <- all.vars(formula)
+ if (length(formula) == 2L) {
+ formula[[3L]] <- formula[[2L]]
+ formula[[2L]] <- 0
+ }
+ form2 <- formula
+ form2[[2L]] <- 0
+ varNamesRHS <- all.vars(form2)
+ mWeights <- missing(weights)
+ pnames <- if (missing(start)) {
+ if (!is.null(attr(data, "parameters"))) {
+ names(attr(data, "parameters"))
+ }
+ else {
+ cll <- formula[[length(formula)]]
+ func <- get(as.character(cll[[1L]]))
+ if (!is.null(pn <- attr(func, "pnames")))
+ as.character(as.list(match.call(func, call = cll))[-1L][pn])
+ }
+ }
+ else names(start)
+ env <- environment(formula)
+ if (is.null(env))
+ env <- parent.frame()
+ if (length(pnames))
+ varNames <- varNames[is.na(match(varNames, pnames))]
+ lenVar <- function(var) tryCatch(length(eval(as.name(var),
+ data, env)), error = function(e) -1)
+ if (length(varNames)) {
+ n <- sapply(varNames, lenVar)
+ if (any(not.there <- n == -1)) {
+ nnn <- names(n[not.there])
+ if (missing(start)) {
+ if (algorithm == "plinear")
+ stop("no starting values specified")
+ warning("No starting values specified for some parameters.\n",
+ "Initializing ", paste(sQuote(nnn), collapse = ", "),
+ " to '1.'.\n", "Consider specifying 'start' or using a selfStart model",
+ domain = NA)
+ start <- setNames(as.list(rep(1, length(nnn))),
+ nnn)
+ varNames <- varNames[i <- is.na(match(varNames,
+ nnn))]
+ n <- n[i]
+ }
+ else stop(gettextf("parameters without starting value in 'data': %s",
+ paste(nnn, collapse = ", ")), domain = NA)
+ }
+ }
+ else {
+ if (length(pnames) && any((np <- sapply(pnames, lenVar)) ==
+ -1)) {
+ message(sprintf(ngettext(sum(np == -1), "fitting parameter %s without any variables",
+ "fitting parameters %s without any variables"),
+ paste(sQuote(pnames[np == -1]), collapse = ", ")),
+ domain = NA)
+ n <- integer()
+ }
+ else stop("no parameters to fit")
+ }
+ respLength <- length(eval(formula[[2L]], data, env))
+ if (length(n) > 0L) {
+ varIndex <- n%%respLength == 0
+ if (is.list(data) && diff(range(n[names(n) %in% names(data)])) >
+ 0) {
+ mf <- data
+ if (!missing(subset))
+ warning("argument 'subset' will be ignored")
+ if (!missing(na.action))
+ warning("argument 'na.action' will be ignored")
+ if (missing(start))
+ start <- getInitial(formula, mf)
+ startEnv <- new.env(hash = FALSE, parent = environment(formula))
+ for (i in names(start)) assign(i, start[[i]], envir = startEnv)
+ rhs <- eval(formula[[3L]], data, startEnv)
+ n <- NROW(rhs)
+ wts <- if (mWeights)
+ rep(1, n)
+ else eval(substitute(weights), data, environment(formula))
+ }
+ else {
+ mf$formula <- as.formula(paste("~", paste(varNames[varIndex],
+ collapse = "+")), env = environment(formula))
+ mf$start <- mf$control <- mf$algorithm <- mf$trace <- mf$model <- NULL
+ mf$lower <- mf$upper <- NULL
+ mf[[1L]] <- as.name("model.frame")
+ mf <- eval.parent(mf)
+ n <- nrow(mf)
+ mf <- as.list(mf)
+ wts <- if (!mWeights)
+ model.weights(mf)
+ else rep(1, n)
+ }
+ if (any(wts < 0 | is.na(wts)))
+ stop("missing or negative weights not allowed")
+ }
+ else {
+ varIndex <- logical()
+ mf <- list(0)
+ wts <- numeric()
+ }
+ if (missing(start))
+ start <- getInitial(formula, mf)
+ for (var in varNames[!varIndex]) mf[[var]] <- eval(as.name(var),
+ data, env)
+ varNamesRHS <- varNamesRHS[varNamesRHS %in% varNames[varIndex]]
+ m <- switch(algorithm, plinear = nlsModel.plinear(formula,
+ mf, start, wts), port = nlsModel(formula, mf, start,
+ wts, upper), nlsModel(formula, mf, start, wts))
+ ctrl <- nls.control()
+ if (!missing(control)) {
+ control <- as.list(control)
+ ctrl[names(control)] <- control
+ }
+ if (algorithm != "port") {
+ if (!missing(lower) || !missing(upper))
+ warning("upper and lower bounds ignored unless algorithm = \"port\"")
+ convInfo <- .Call(C_nls_iter, m, ctrl, trace)
+ nls.out <- list(m = m, convInfo = convInfo, data = substitute(data),
+ call = match.call())
+ }
+ else {
+ pfit <- nls_port_fit(m, start, lower, upper, control,
+ trace, give.v = TRUE)
+ iv <- pfit[["iv"]]
+ msg.nls <- port_msg(iv[1L])
+ conv <- (iv[1L] %in% 3:6)
+ if (!conv) {
+ msg <- paste("Convergence failure:", msg.nls)
+ if (ctrl$warnOnly)
+ warning(msg)
+ else stop(msg)
+ }
+ v. <- port_get_named_v(pfit[["v"]])
+ cInfo <- list(isConv = conv, finIter = iv[31L], finTol = v.[["NREDUC"]],
+ nEval = c(`function` = iv[6L], gradient = iv[30L]),
+ stopCode = iv[1L], stopMessage = msg.nls)
+ cl <- match.call()
+ cl$lower <- lower
+ cl$upper <- upper
+ nls.out <- list(m = m, data = substitute(data), call = cl,
+ convInfo = cInfo, convergence = as.integer(!conv),
+ message = msg.nls)
+ }
+ nls.out$call$algorithm <- algorithm
+ nls.out$call$control <- ctrl
+ nls.out$call$trace <- trace
+ nls.out$na.action <- attr(mf, "na.action")
+ nls.out$dataClasses <- attr(attr(mf, "terms"), "dataClasses")[varNamesRHS]
+ if (model)
+ nls.out$model <- mf
+ if (!mWeights)
+ nls.out$weights <- wts
+ nls.out$control <- control
+ class(nls.out) <- "nls"
+ nls.out
+ if(is.null(vcov)) {
+ se <- vcov(nls.out)
+ } else {
+ if (is.function(vcov))
+ se <- vcov(nls.out)
+ else
+ se <- vcov
+ }
+ nls.out = list(nls.out,vHaC = se)
+ nls.out
+
+}
\ No newline at end of file
Added: pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Literature/Thumbs.db
===================================================================
(Binary files differ)
Property changes on: pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Literature/Thumbs.db
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/lmi.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/lmi.R (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/lmi.R 2013-08-03 10:10:40 UTC (rev 2713)
@@ -0,0 +1,76 @@
+lmi <- function (formula, data,vcov = NULL, subset, weights, na.action, method = "qr",
+ model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE,
+ contrasts = NULL, offset, ...)
+{
+ ret.x <- x
+ ret.y <- y
+ cl <- match.call()
+ mf <- match.call(expand.dots = FALSE)
+ m <- match(c("formula", "data", "subset", "weights", "na.action",
+ "offset"), names(mf), 0L)
+ mf <- mf[c(1L, m)]
+ mf$drop.unused.levels <- TRUE
+ mf[[1L]] <- as.name("model.frame")
+ mf <- eval(mf, parent.frame())
+ if (method == "model.frame")
+ return(mf)
+ else if (method != "qr")
+ warning(gettextf("method = '%s' is not supported. Using 'qr'",
+ method), domain = NA)
+ mt <- attr(mf, "terms")
+ y <- model.response(mf, "numeric")
+ w <- as.vector(model.weights(mf))
+ if (!is.null(w) && !is.numeric(w))
+ stop("'weights' must be a numeric vector")
+ offset <- as.vector(model.offset(mf))
+ if (!is.null(offset)) {
+ if (length(offset) != NROW(y))
+ stop(gettextf("number of offsets is %d, should equal %d (number of observations)",
+ length(offset), NROW(y)), domain = NA)
+ }
+ if (is.empty.model(mt)) {
+ x <- NULL
+ z <- list(coefficients = if (is.matrix(y)) matrix(, 0,
+ 3) else numeric(), residuals = y, fitted.values = 0 *
+ y, weights = w, rank = 0L, df.residual = if (!is.null(w)) sum(w !=
+ 0) else if (is.matrix(y)) nrow(y) else length(y))
+ if (!is.null(offset)) {
+ z$fitted.values <- offset
+ z$residuals <- y - offset
+ }
+ }
+ else {
+ x <- model.matrix(mt, mf, contrasts)
+ z <- if (is.null(w))
+ lm.fit(x, y, offset = offset, singular.ok = singular.ok,
+ ...)
+ else lm.wfit(x, y, w, offset = offset, singular.ok = singular.ok,
+ ...)
+ }
+ class(z) <- c(if (is.matrix(y)) "mlm", "lm")
+ z$na.action <- attr(mf, "na.action")
+ z$offset <- offset
+ z$contrasts <- attr(x, "contrasts")
+ z$xlevels <- .getXlevels(mt, mf)
+ z$call <- cl
+ z$terms <- mt
+ if (model)
+ z$model <- mf
+ if (ret.x)
+ z$x <- x
+ if (ret.y)
+ z$y <- y
+ if (!qr)
+ z$qr <- NULL
+ #z
+ if(is.null(vcov)) {
+ se <- vcov(z)
+ } else {
+ if (is.function(vcov))
+ se <- vcov(z)
+ else
+ se <- vcov
+ }
+ z = list(z,vHaC = se)
+ z
+}
More information about the Returnanalytics-commits
mailing list