diff options
Diffstat (limited to 'R/poly.R')
| -rw-r--r-- | R/poly.R | 196 |
1 files changed, 196 insertions, 0 deletions
diff --git a/R/poly.R b/R/poly.R new file mode 100644 index 0000000..528a599 --- /dev/null +++ b/R/poly.R @@ -0,0 +1,196 @@ +#' Polynomial model with identifiable parameters +#' +#' Creates a polynomial model with identifiable coefficients +#' +#' @param A autoregressive coefficients +#' @param B,F1 coefficients of the numerator and denominator respectively +#' of the deterministic model between the input and output +#' @param C,D coefficients of the numerator and denominator respectively +#' of the stochastic model +#' @param ioDelay the delay in the input-output channel +#' @param Ts sampling interval +#' @param noiseVar variance of the white noise source (Default=\code{1}) +#' @param intNoise Logical variable indicating presence or absence of integrator +#' in the noise channel (Default=\code{FALSE}) +#' @param unit time unit (Default=\code{"seconds"}) +#' +#' @details +#' Discrete-time polynomials are of the form +#' \deqn{ +#' A(q^{-1}) y[k] = \frac{B(q^{-1})}{F1(q^{-1})} u[k] + +#' \frac{C(q^{-1})}{D(q^{-1})} e[k] +#' } +#' +#' @examples +#' # define output-error model +#' mod_oe <- idpoly(B=c(0.6,-0.2),F1=c(1,-0.5),ioDelay = 2,Ts=0.1, +#' noiseVar = 0.1) +#' +#' # define box-jenkins model with unit variance +#' B <- c(0.6,-0.2) +#' C <- c(1,-0.3) +#' D <- c(1,1.5,0.7) +#' F1 <- c(1,-0.5) +#' mod_bj <- idpoly(1,B,C,D,F1,ioDelay=1) +#' +#' @export +idpoly <- function(A=1,B=1,C=1,D=1,F1=1,ioDelay=0,Ts=1, + noiseVar=1,intNoise = F,unit = c("seconds","minutes", + "hours","days")[1]){ + Bindex <- which.max(B!=0) + ioDelay <- ifelse(Bindex-1==0,ioDelay,Bindex-1) + B <- B[Bindex:length(B)] + out <- list(A= A,B=B,C=C,D=D,F1=F1,ioDelay = ioDelay,Ts=Ts, + noiseVar=noiseVar,unit=unit,intNoise = intNoise) + out$type <- typecheck(out) + class(out) <- "idpoly" + return(out) +} + +typecheck <- function(x){ + y <- lapply(x[1:5],checkUnity) + if(y$A){ + out <- if(y$C && y$D) if(y$F1) "fir" else "oe" else "bj" + } else{ + if(y$D && y$F1){ + if(x$intNoise){ + out <- if(y$C) "ari" else "arima" + } else{ + out <- if(y$C) "ar" else "arma" + } + if(!y$B) out <- paste(out,"x",sep="") + } else{ + out <- "polynomial" + } + } + out +} + +checkUnity <- function(x){ + out <- if(length(x)==1 && x==1) TRUE else FALSE +} + +#' @export +print.idpoly <- function(x,se=NULL,dig=3,...){ + main <- paste("Discrete-time",toupper(x$type),"model:") + if(x$type=="oe" || x$type=="bj"){ + main <- paste(main,"y[k] = B(z)/F(z) u[k] +") + if(x$type=="bj"){ + polyExp <- if(!checkUnity(x$C)) "C(z)/D(z)" else "1/D(z)" + if(x$intNoise==T) polyExp <- paste(polyExp,"1/(1-z^{-1})") + main <- paste(main,polyExp,"e[k] \n\n") + } + } else{ + main <- paste(main,"A(z)y[k] =") + if(!checkUnity(x$B)) main <- paste(main,"B(z)u[k] +") + if(checkUnity(x$C)){ + Cexp <- if(!x$intNoise) "" else "1/(1-z^{-1})" + }else{ + Cexp <- "C(z)" + if(x$intNoise) Cexp <- paste(Cexp,"/(1-z^{-1})",sep="") + } + main <- paste(main,Cexp,"e[k] \n\n") + } + cat(main) + + # Printing Standard error sequence + j=1 + print_se <- function(se){ + if(!is.null(se)){ + cat(" (+/- ",round(se[j],dig),") ",sep = "") + j <<- j+1 + } + } + + if(length(x$A)>1){ + cat("A(q^{-1}) = ") + for(i in seq_along(x$A)){ + if(i-1==0){ + cat(round(x$A[i],dig)) + } else{ + if(x$A[i]>0) cat(" + ") else cat("- ") + + if(!(abs(x$A[i])==1)) cat(abs(round(x$A[i],dig))) + print_se(se) + cat("q^{-",i-1,"}",sep="") + } + cat("\t") + } + cat("\n") + } + + cat("B(q^{-1}) = ") + for(i in seq_along(x$B)){ + if(i+x$ioDelay-1==0){ + cat(round(x$B[i],dig)) + } else{ + + if(!((x$ioDelay!=0) && (i==1))){ + if(x$B[i]>0) cat(" + ") else cat("- ") + } else{ + if(x$B[i]<0) cat("-") + } + + if(!(abs(x$B[i])==1)) cat(abs(round(x$B[i],dig))) + print_se(se) + cat("q^{-",i+x$ioDelay-1,"}",sep="") + } + cat("\t") + } + cat("\n") + + if(length(x$C)>1){ + cat("C(q^{-1}) = ") + for(i in seq_along(x$C)){ + if(i-1==0){ + cat(round(x$C[i],dig)) + } else{ + if(x$C[i]>0) cat(" + ") else cat("- ") + + if(!(abs(x$C[i])==1)) cat(abs(round(x$C[i],dig))) + print_se(se) + cat("q^{-",i-1,"}",sep="") + } + cat("\t") + } + cat("\n") + } + + if(length(x$D)>1){ + cat("D(q^{-1}) = ") + for(i in seq_along(x$D)){ + if(i-1==0){ + cat(round(x$D[i],dig)) + } else{ + if(x$D[i]>0) cat(" + ") else cat("- ") + + if(!(abs(x$D[i])==1)) cat(abs(round(x$D[i],dig))) + print_se(se) + cat("q^{-",i-1,"}",sep="") + } + cat("\t") + } + cat("\n") + } + + if(length(x$F1)>1){ + cat("F(q^{-1}) = ") + for(i in seq_along(x$F1)){ + if(i-1==0){ + cat(round(x$F1[i],dig)) + } else{ + if(x$F1[i]>0) cat(" + ") else cat("- ") + + if(!(abs(x$F1[i])==1)) cat(abs(round(x$F1[i],dig))) + print_se(se) + cat("q^{-",i-1,"}",sep="") + } + cat("\t") + } + } + cat("\n") +} + +params <- function(x){ + c(x$A[-1],x$B,x$C[-1],x$D[-1],x$F1[-1]) +}
\ No newline at end of file |