diff options
| -rw-r--r-- | R/estpoly.R | 39 |
1 files changed, 28 insertions, 11 deletions
diff --git a/R/estpoly.R b/R/estpoly.R index 7664264..3074b4c 100644 --- a/R/estpoly.R +++ b/R/estpoly.R @@ -282,7 +282,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^(-5); sumSqRatio <- 1000; lambda <- 1 + tol <- 10^(-5); sumSqRatio <- 1000; lambda <- 10 eout <- matrix(rep(0,N+2*n)) reg <- function(i) { @@ -292,30 +292,47 @@ armax <- function(x,order=c(0,1,1,0)){ # Initialize Algorithm i = 0 - theta <- matrix(runif(na+nb+nc,min=-0.5,max=0.5)) + 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%*%theta + e <- Y-X%*%theta0 + sumsq0 <- sum(e^2) + updateJ <- 1 while (sumSqRatio > tol){ - sumsq0 <- sum(e^2) - # Compute gradient - eout <- matrix(c(rep(0,n),e[,])) - X <- t(sapply(n+1:(N+n),reg)) - filt1 <- Arma(b=1,a=c(1,theta[(na+nb+1:nc)])) - grad <- apply(X,2,filter,filt=filt1) + + 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 <- theta + 1/N*Hinv%*%t(grad)%*%e + 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) + # print(sumsq);print(sumSqRatio) i=i+1 } # print(sumSqRatio) |