diff options
Diffstat (limited to 'R')
-rw-r--r-- | R/estpoly.R | 25 |
1 files changed, 15 insertions, 10 deletions
diff --git a/R/estpoly.R b/R/estpoly.R index 83aa204..f0382fd 100644 --- a/R/estpoly.R +++ b/R/estpoly.R @@ -177,6 +177,7 @@ arx <- function(x,order=c(0,1,0)){ residuals=(Y-X%*%coef)[1:N,],call=match.call(),input=u) } +#' @export armax <- function(x,order=c(0,1,1,0)){ require(signal) require(MASS) @@ -190,7 +191,7 @@ 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^(-3); sumsq <- 10^3; i = 0 + tol <- 10^(-3); sumSqRatio <- 1000 eout <- matrix(rep(0,N+2*n)) reg <- function(i) { @@ -198,15 +199,15 @@ armax <- function(x,order=c(0,1,1,0)){ matrix(c(-yout[i-1:na,],uout[v,],eout[i-1:nc,])) } - while ((sumsq > tol) && (i<20)){ - if(i==0){ - theta <- matrix(runif(na+nb+nc,min=-0.2,max=0.2)) - } - - # Generate Residuals from previous params - X <- t(sapply(n+1:(N+n),reg)) - Y <- yout[n+1:(N+n),,drop=F] - e <- Y-X%*%theta + # Initialize Algorithm + i = 0 + theta <- matrix(runif(na+nb+nc,min=-0.2,max=0.2)) + X <- t(sapply(n+1:(N+n),reg)) + Y <- yout[n+1:(N+n),,drop=F] + e <- Y-X%*%theta + + while ((sumSqRatio > tol) && (i<50)){ + sumsq0 <- sum(e^2) # Compute gradient eout <- matrix(c(rep(0,n),e[,])) @@ -217,7 +218,11 @@ armax <- function(x,order=c(0,1,1,0)){ # Update Parameters theta <- theta + ginv(grad)%*%e + # Update residuals + e <- Y-X%*%theta sumsq <- sum(e^2) + + sumSqRatio <- abs(sumsq0-sumsq)/sumsq0 i=i+1 } model <- idpoly(A = c(1,theta[1:na]),B = theta[na+1:nb], |