R/ftf.fit.R

In spa: Implements The Sequential Predictions Algorithm

"ftf.fit" <-function(y,S,type="soft",global=mean(y),control,GCV=FALSE,F1=NULL){
  val=dim(as.matrix(y))
  m=val[1]
  n<-dim(S)[1]
  if(n==m){
    S=S/apply(S,1,sum)
    prn=as.vector(S%*%y)
    if(GCV)
      return(list(prn,sum(diag(S))))
    return(prn)
  }
  prn<-rep(global,n)
  SUL=as.matrix(S[(m+1):n,1:m])
  SUU=as.matrix(S[(m+1):n,(m+1):n])
  if(dim(SUL)[2]==1){
    a1=sum(SUL)
    if(a1==0){
      S=S/apply(S,1,sum)
      vec=rep(0,n)
      vec[1:m]=y
      prn=as.vector(S%*%vec)
      prn[n]=global
      if(GCV)
        return(list(prn,sum(diag(S[1:m,1:m]))))
      return(prn)
    }
    a=1/(1-SUU[1,1])* sum(SUL*y)
    b=as.vector(S[1:m,1:m]%*%y+S[1:m,(m+1):n]*a)
    prn=c(b,a)
    if(GCV)
      return(list(prn,sum(diag(S[1:m,1:m]))))
    return(prn)
  }
  ignore<-rep(FALSE,n)
  v=apply(SUL,1,sum)
  ignore[(m+1):n]<-v==0
  if(sum(ignore)>0){
    n<-sum(!ignore)
    S=S[!ignore,!ignore]
    if(m==n){
      S=S/apply(S,1,sum)
      prn[!ignore]=as.vector(S%*%y)
      if(GCV)
        return(list(prn,sum(diag(S[1:m,1:m]))))
      return(prn)
    }
    SUL=as.matrix(S[(m+1):n,1:m])
    SUU=as.matrix(S[(m+1):n,(m+1):n])
    if(dim(SUL)[2]==1){
      a1=sum(SUL)
      if(a1==0){
        S=S/apply(S,1,sum)
        vec=rep(0,n)
        vec[1:m]=y
        prn=as.vector(S%*%vec)
        prn[n]=global
        if(GCV)
          return(list(prn,sum(diag(S[1:m,1:m]))))
        return(prn)
      }
      a=1/(1-SUU[1,1])* sum(SUL*y)
      b=as.vector(S[1:m,1:m]%*%y+S[1:m,(m+1):n]*a)
      prn[!ignore]=c(b,a)
      if(GCV){
        V=S[1:m,(m+1):n]%*%SUL
        V=V*1/(1-SUU[1,1]) 
        A=S[1:m,1:m]+V
        return(list(prn,sum(diag(A))))
      }
      return(prn)
    }
  }
  M=try(solve(diag(n-m)-SUU,SUL),TRUE)
  if(class(M)=="try-error"){
    if(control$warn)
      warning("Matrix Inverse Problem ... Using ginv instead of solve\n")
    M=ginv(diag(n-m)-SUU)%*%SUL
  }
  a=as.vector(M%*%y)
  A=S[1:m,1:m]+S[1:m,(m+1):n]%*%M
  b=as.vector(A%*%y)
  if(type=="soft"){
    prn[!ignore]=c(b,a)
    if(GCV)
      return(list(prn,sum(diag(A))))
    return(prn)
  }
  if(is.null(F1)){
    pr1=c(b,a)
  }else{
    F1=F1[,1:m]/apply(F1[1:m,],2,sum)
    pr1=F1%*%y
  }
  pr1=pr1[-c(1:m)]
  eps=control$eps
  maxiter=control$maxiter
  for(i in 1:maxiter){
    pro=pr1
    h=c(y,round(pr1))
    h=as.vector(S%*%h)
    pr1=h[-c(1:m)]
    tol=sum( (pr1-pro)^2)
    if(tol<eps)
      break
  }
  prn[!ignore]=h
  if(GCV){
    return(list(prn,sum(diag(A))))
  }
  return(prn)
}