diff options
| -rw-r--r-- | NAMESPACE | 9 | ||||
| -rw-r--r-- | R/estpoly.R | 101 | ||||
| -rw-r--r-- | man/arx.Rd (renamed from man/estARX.Rd) | 8 |
3 files changed, 63 insertions, 55 deletions
@@ -17,21 +17,22 @@ S3method(plot,idframe) S3method(plot,idfrd) S3method(plot,impulseest) S3method(predict,detrend) -S3method(predict,estARX) +S3method(predict,estPoly) S3method(print,idpoly) -S3method(print,summary.estARX) +S3method(print,summary.estPoly) S3method(print,summary.idframe) S3method(print,tf) S3method(sim,default) S3method(sim,idpoly) -S3method(summary,estARX) +S3method(summary,estPoly) S3method(summary,idframe) S3method(time,idframe) export("inputNames<-") export("outputNames<-") +export(arx) export(dataSlice) export(detrend) -export(estARX) +export(estPoly) export(etfe) export(idframe) export(idfrd) diff --git a/R/estpoly.R b/R/estpoly.R index 79b6f5c..66176dc 100644 --- a/R/estpoly.R +++ b/R/estpoly.R @@ -1,4 +1,55 @@ #' @export +estPoly <- function(coefficients,vcov,sigma,df,fitted.values, + residuals,call,time,input){ + out <- list(coefficients= coefficients,vcpv= vcov,sigma = sigma, + df= df,fitted.values=fitted.values, + residuals= residuals,call= call,time=time,input=input) +} + +#' @export +summary.estPoly <- function(object) +{ + coefs <- c(coef(object)$A[-1],coef(object)$B) + se <- sqrt(diag(object$vcov)) + tval <- coefs / se + TAB <- cbind(Estimate = coefs, + StdErr = se, + t.value = tval, + p.value = 2*pt(-abs(tval), df=object$df)) + na <- length(coef(object)$A) - 1; nk <- coef(object)$ioDelay; + nb <- length(coef(object)$B) - nk + + rownames(TAB) <- rep("a",nrow(TAB)) + for(i in 1:na) rownames(TAB)[i] <- paste("a",i,sep="") + for(j in (na+1):nrow(TAB)) { + rownames(TAB)[j] <- paste("b",j-na-1+nk,sep="") + } + res <- list(call=object$call,coefficients=TAB,sigma=object$sigma, + df=object$df) + class(res) <- "summary.estARX" + res +} + +#' @export +print.summary.estPoly <- function(object){ + cat("Discrete-time ARX model: A(q^{-1})y[k] = B(q^{-1})u[k] + e[k] \n") + cat("Call: ");print(object$call);cat("\n\n") + + print(coef(object)) + cat(paste("\nsigma:",format(object$sigma,digits=4))) + cat(paste("\nDoF:",object$df)) +} + +#' @export +predict.estPoly <- function(model,newData=NULL){ + if(is.null(newdata)){ + return(fitted(model)) + } else{ + return(sim(coef(model),newdata$input)) + } +} + +#' @export plot.estPoly <- function(model,newdata=NULL){ require(ggplot2) @@ -84,12 +135,12 @@ residplot <- function(model,newdata=NULL){ #' #' @examples #' data(arxsim) -#' model <- estARX(data,c(2,1,1)) +#' model <- arx(data,c(2,1,1)) #' summary(model) # obtain estimates and their covariances #' plot(model) # plot the predicted and actual responses #' #' @export -estARX <- function(data,order=c(0,1,0)){ +arx <- function(data,order=c(0,1,0)){ y <- as.matrix(data$output) u <- as.matrix(data$input); N <- dim(y)[1] na <- order[1];nb <- order[2]; nk <- order[3] @@ -114,53 +165,9 @@ estARX <- function(data,order=c(0,1,0)){ model <- arx(A = c(1,coef[1:na]),B = coef[na+1:(nb+1)],ioDelay = nk) time <- seq(from=data$t.start,to=data$t.end,by=data$Ts) - est <- list(coefficients = model,vcov = vcov, sigma = sqrt(sigma2), + est <- estPoly(coefficients = model,vcov = vcov, sigma = sqrt(sigma2), df = df,fitted.values=(X%*%coef)[1:N,], residuals=(Y-X%*%coef)[1:N,],call=match.call(), time=time,input=u) - class(est) <- c("estARX","estPoly") est -} - -#' @export -predict.estARX <- function(model,newdata=NULL){ - if(is.null(newdata)){ - return(fitted(model)) - } else{ - return(sim(coef(model),newdata$input)) - } -} - -#' @export -summary.estARX <- function(object) -{ - coefs <- c(coef(object)$A[-1],coef(object)$B) - se <- sqrt(diag(object$vcov)) - tval <- coefs / se - TAB <- cbind(Estimate = coefs, - StdErr = se, - t.value = tval, - p.value = 2*pt(-abs(tval), df=object$df)) - na <- length(coef(object)$A) - 1; nk <- coef(object)$ioDelay; - nb <- length(coef(object)$B) - nk - - rownames(TAB) <- rep("a",nrow(TAB)) - for(i in 1:na) rownames(TAB)[i] <- paste("a",i,sep="") - for(j in (na+1):nrow(TAB)) { - rownames(TAB)[j] <- paste("b",j-na-1+nk,sep="") - } - res <- list(call=object$call,coefficients=TAB,sigma=object$sigma, - df=object$df) - class(res) <- "summary.estARX" - res -} - -#' @export -print.summary.estARX <- function(object){ - cat("Discrete-time ARX model: A(q^{-1})y[k] = B(q^{-1})u[k] + e[k] \n") - cat("Call: ");print(object$call);cat("\n\n") - - print(coef(object)) - cat(paste("\nsigma:",format(object$sigma,digits=4))) - cat(paste("\nDoF:",object$df)) }
\ No newline at end of file diff --git a/man/estARX.Rd b/man/arx.Rd index 34c0795..62b9cd1 100644 --- a/man/estARX.Rd +++ b/man/arx.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2 (4.1.1): do not edit by hand % Please edit documentation in R/estpoly.R -\name{estARX} -\alias{estARX} +\name{arx} +\alias{arx} \title{Estimate ARX Models} \usage{ -estARX(data, order = c(0, 1, 0)) +arx(data, order = c(0, 1, 0)) } \arguments{ \item{data}{an object of class \code{idframe}} @@ -48,7 +48,7 @@ using the \code{\link{detrend}} function. } \examples{ data(arxsim) -model <- estARX(data,c(2,1,1)) +model <- arx(data,c(2,1,1)) summary(model) # obtain estimates and their covariances plot(model) # plot the predicted and actual responses } |