modifying estimation routines to desired format

6 files changed, 74 insertions, 95 deletions

diff --git a/NAMESPACE b/NAMESPACE
index 2764a50..20f4528 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -12,20 +12,20 @@ S3method(outputData,default)
 S3method(outputData,idframe)
 S3method(outputNames,default)
 S3method(outputNames,idframe)
-S3method(plot,estPoly)
+S3method(plot,estpoly)
 S3method(plot,idframe)
 S3method(plot,idfrd)
 S3method(plot,impulseest)
 S3method(predict,detrend)
-S3method(predict,estPoly)
-S3method(print,estPoly)
+S3method(predict,estpoly)
+S3method(print,estpoly)
 S3method(print,idpoly)
-S3method(print,summary.estPoly)
+S3method(print,summary.estpoly)
 S3method(print,summary.idframe)
 S3method(print,tf)
 S3method(sim,default)
 S3method(sim,idpoly)
-S3method(summary,estPoly)
+S3method(summary,estpoly)
 S3method(summary,idframe)
 S3method(time,idframe)
 export("inputNames<-")
@@ -34,8 +34,9 @@ export(armax)
 export(arx)
 export(dataSlice)
 export(detrend)
-export(estPoly)
+export(estpoly)
 export(etfe)
+export(getcov)
 export(idframe)
 export(idfrd)
 export(idinput)
@@ -47,6 +48,7 @@ export(misdata)
 export(nInputSeries)
 export(nOutputSeries)
 export(oe)
+export(optimOptions)
 export(outputData)
 export(outputNames)
 export(read.idframe)
diff --git a/R/estUtil.R b/R/estUtil.R
index ba9f6d1..d17ccbf 100644
--- a/R/estUtil.R
+++ b/R/estUtil.R
@@ -1,18 +1,46 @@
+armaxGrad <- function(theta,e,dots){
+  y <- dots[[1]]; u <- dots[[2]]; order <- dots[[3]]; N <- dots[[4]]
+  na <- order[1];nb <- order[2]; nc <- order[3]; nk <- order[4]
+  nb1 <- nb+nk-1 ; n <- max(na,nb1,nc)
+  
+  if(is.null(e)){
+    eout <- matrix(rep(0,N+2*n))
+  } else{
+    eout <- matrix(c(rep(0,n),e[,]))
+  }
+  
+  reg <- function(i) {
+    if(nk==0) v <- i-0:(nb-1) else v <- i-nk:nb1
+    matrix(c(-y[i-1:na,],u[v,],eout[i-1:nc,]))
+  }
+  
+  X <- t(sapply(n+1:(N+n),reg))
+  Y <- y[n+1:(N+n),,drop=F]
+  l <- list(X=X,Y=Y)
+  
+  if(!is.null(e)){
+    filt1 <- Arma(b=1,a=c(1,theta[(na+nb+1:nc)]))
+    grad <- apply(X,2,filter,filt=filt1) 
+    l$grad <- grad
+  }
+  
+  return(l)
+}
+
 # Implementation of the Levenberg Marquardt Algorithm
-levbmqdt <- function(y,u,n,reg,theta0,opt=optimOptions()){
+levbmqdt <- function(...,obj,theta0,N,opt=optimOptions()){
+  dots <- list(...)
   
   # Optimization Parameters
   tol <- opt$tol; maxIter <- opt$maxIter
   d <- opt$adv$LMinit; mu <- opt$adv$LMstep
   
-  N <- dim(y,1); df <- N - length(theta0)
+  df <- N - dim(theta0)[1]
   
   # Initialize Algorithm
   i <- 0
-  eout <- matrix(rep(0,N+2*n))
-  X <- t(sapply(n+1:(N+n),reg,y,u,e))
-  Y <- yout[n+1:(N+n),,drop=F]
-  e <- Y-X%*%theta0
+  l <- obj(theta=theta0,e=NULL,dots)
+  e <- l$Y-l$X%*%theta0
   sumsq0 <- sum(e^2)
   # parameter indicating whether to update gradient in the iteration
   update <- 1 
@@ -24,39 +52,35 @@ levbmqdt <- function(y,u,n,reg,theta0,opt=optimOptions()){
     if(update ==1){
       countObj <- countObj+1
       # Update gradient
-      eout <- matrix(c(rep(0,n),e[,]))
-      X <- t(sapply(n+1:(N+n),reg,y,u,e))
-      filt1 <- Arma(b=1,a=c(1,theta0[(na+nb+1:nc)]))
-      grad <- apply(X,2,filter,filt=filt1) 
+      l <- obj(theta0,e,dots)
     }
     
     # Update Parameters
-    H <- 1/N*(t(grad)%*%grad) + d*diag(na+nb+nc)
+    H <- 1/N*(t(l$grad)%*%l$grad) + d*diag(dim(theta0)[1])
     Hinv <- solve(H)
-    theta <- theta0 + 1/N*Hinv%*%t(grad)%*%e
+    theta <- theta0 + 1/N*Hinv%*%t(l$grad)%*%e
     
     # Update residuals
-    e <- Y-X%*%theta
+    e <- l$Y-l$X%*%theta
     sumsq <- sum(e^2)
+    sumSqRatio <- (sumsq0-sumsq)/sumsq0*100
     
     # If sum square error with the updated parameters is less than the 
     # previous one, the updated parameters become the current parameters
     # and the damping coefficient is reduced by a factor of mu
-    if(sumsq<sumsq0){
+    if(sumSqRatio > 0){
       d <- d/mu
       theta0 <- theta
       sumsq0 <- sumsq
       update <- 1
     } else{ # increase damping coefficient by a factor of mu
-      d <- d*v
+      d <- d*mu
       update <- 0
     }
     
-    sumSqRatio <- abs(sumsq0-sumsq)/sumsq0*100
-    #     print(sumsq);print(sumSqRatio)
-    
-    if(sumSqRatio < tol || i > maxIter)
+    if((sumSqRatio < tol) || (i == maxIter)){
       break
+    }
   }
   
   if(sumSqRatio < tol){
@@ -74,7 +98,12 @@ levbmqdt <- function(y,u,n,reg,theta0,opt=optimOptions()){
 }
 
 #' @export
-optimOptions <- function(tol=0.01,maxIter=50,LMinit=0.1,LMstep=2){
+optimOptions <- function(tol=0.01,maxIter=20,LMinit=10,LMstep=2){
   return(list(tol=tol,maxIter= maxIter, adv= list(LMinit=LMinit,
                                                   LMstep=LMstep)))
+}
+
+#' @export
+getcov <- function(sys){
+  sys$stats$vcov
 }
\ No newline at end of file
diff --git a/R/estpoly.R b/R/estpoly.R
index 63b5597..679281d 100644
--- a/R/estpoly.R
+++ b/R/estpoly.R
@@ -1,13 +1,14 @@
 #' @export
-estpoly <- function(coefficients,vcov,sigma,df,fitted.values,
-                    residuals,call,input){
-  out <- list(coefficients= coefficients,vcov= vcov,sigma = sigma,
-              df= df,fitted.values=fitted.values, 
-              residuals= residuals,call= call,input=input)
+estpoly <- function(model,fitted.values,residuals,options=NULL,
+                    call,stats,termination=NULL,datainfo){
+  out <- list(model=model,fitted.values=fitted.values,residuals=residuals,
+              datainfo=datainfo,call=call,stats=stats,options=options,
+              termination=termination)
   class(out) <- "estpoly"
   out
 }
 
+
 #' @export
 print.estpoly <- function(est,...){
   print(summary(est),...)
@@ -16,7 +17,7 @@ print.estpoly <- function(est,...){
 #' @export
 summary.estpoly <- function(object)
 {
-  model <- coef(object)
+  model <- estpoly$sys
   if(model$type=="arx"||model$type=="armax"){
     coefs <- c(model$A[-1],model$B)
     na <- length(model$A) - 1; nk <- model$ioDelay; 
@@ -31,12 +32,7 @@ summary.estpoly <- function(object)
     nb <- length(model$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))
+  se <- sqrt(diag(getcov(object)))
   
   rownames(TAB) <- rep("a",nrow(TAB))
   
@@ -282,69 +278,21 @@ armax <- function(x,order=c(0,1,1,0)){
   padZeros <- function(x,n) c(rep(0,n),x,rep(0,n))
   yout <- apply(y,2,padZeros,n=n)
   uout <- apply(u,2,padZeros,n=n)
-  tol <- 10^(-5); sumSqRatio <- 1000; lambda <- 10
-  eout <- matrix(rep(0,N+2*n))
   
-  reg <- function(i) {
+  reg <- function(i,y,u,e) {
     if(nk==0) v <- i-0:(nb-1) else v <- i-nk:nb1
-    matrix(c(-yout[i-1:na,],uout[v,],eout[i-1:nc,]))
+    matrix(c(-y[i-1:na,],u[v,],e[i-1:nc,]))
   }
-  
-  # Initialize Algorithm
-  i = 0
   theta0 <- matrix(rnorm(na+nb+nc)) # current parameters
-  X <- t(sapply(n+1:(N+n),reg))
-  Y <- yout[n+1:(N+n),,drop=F]
-  e <- Y-X%*%theta0
-  sumsq0 <- sum(e^2)
-  updateJ <- 1
   
-  while (sumSqRatio > tol){
-    
-    if(updateJ ==1){
-      # Compute gradient
-      eout <- matrix(c(rep(0,n),e[,]))
-      X <- t(sapply(n+1:(N+n),reg))
-      filt1 <- Arma(b=1,a=c(1,theta0[(na+nb+1:nc)]))
-      grad <- apply(X,2,filter,filt=filt1) 
-    }
-    
-    # Update Parameters
-    H <- 1/N*(t(grad)%*%grad) + lambda*diag(na+nb+nc)
-    Hinv <- solve(H)
-    theta <- theta0 + 1/N*Hinv%*%t(grad)%*%e
-    
-    # Update residuals
-    e <- Y-X%*%theta
-    sumsq <- sum(e^2)
-    
-    # If sum square error with the updated parameters is less than the 
-    # previous one, the updated parameters become the current parameters
-    # and the damping factor is reduced by a factor of 10
-    if(sumsq<sumsq0){
-      lambda <- lambda/10
-      theta0 <- theta
-      sumsq0 <- sumsq
-      updateJ <- 1
-    } else{ # increase damping factor
-      lambda <- lambda*10
-      updateJ <- 0
-    }
-    
-    sumSqRatio <- abs(sumsq0-sumsq)/sumsq0
-    #     print(sumsq);print(sumSqRatio)
-    i=i+1
-  }
-  # print(sumSqRatio)
-  e <- e[1:N,]
-  sigma2 <- sum(e^2)/df
-  vcov <- sigma2 * Hinv
+  l <- levbmqdt(yout,uout,order,N,obj=armaxGrad,theta0=theta0,N=N)
+  theta <- l$theta
   
   model <- idpoly(A = c(1,theta[1:na]),B = theta[na+1:nb],
                   C = c(1,theta[na+nb+1:nc]),ioDelay = nk,Ts=deltat(x))
   
-  estpoly(coefficients = model,vcov = vcov, sigma = sqrt(sigma2),
-          df = df,fitted.values=y-e, residuals=e,call=match.call(),
+  estpoly(coefficients = model,vcov = l$vcov, sigma = l$sigma,df = df,
+          fitted.values=y, residuals=l$e,call=match.call(),
           input=u)
 }
 
diff --git a/man/armax.Rd b/man/armax.Rd
index d4e5cf5..1a711e3 100644
--- a/man/armax.Rd
+++ b/man/armax.Rd
@@ -14,7 +14,7 @@ armax(x, order = c(0, 1, 1, 0))
 + 1, order of the polynomial C,and the input-output delay respectively}
 }
 \value{
-An object of class \code{estPoly} containing the following elements:
+An object of class \code{estpoly} containing the following elements:
 
 \tabular{ll}{
    \code{coefficients} \tab an \code{idpoly} object containing the
diff --git a/man/arx.Rd b/man/arx.Rd
index 0ed53f6..f9d38b8 100644
--- a/man/arx.Rd
+++ b/man/arx.Rd
@@ -14,7 +14,7 @@ arx(x, order = c(0, 1, 0))
 the input-output delay}
 }
 \value{
-An object of class \code{estPoly} containing the following elements:
+An object of class \code{estpoly} containing the following elements:
 
 \tabular{ll}{
    \code{coefficients} \tab an \code{idpoly} object containing the
diff --git a/man/oe.Rd b/man/oe.Rd
index a893bb6..78d2de7 100644
--- a/man/oe.Rd
+++ b/man/oe.Rd
@@ -14,7 +14,7 @@ oe(x, order = c(1, 1, 0))
 and the input-output delay respectively}
 }
 \value{
-An object of class \code{estPoly} containing the following elements:
+An object of class \code{estpoly} containing the following elements:
 
 \tabular{ll}{
    \code{coefficients} \tab an \code{idpoly} object containing the