R/simpls.R

In ftsa: Functional Time Series Analysis

simpls <- function(Xtrain, Ytrain, Xtest, ncomp = NULL, weight = FALSE, beta = 0.1)
{
    X = scale(Xtrain, center = TRUE, scale = FALSE)
    if (is.vector(Ytrain)){
       Ytrain = matrix(Ytrain, length(Ytrain), 1)
    }
    Y = scale(Ytrain, center = TRUE, scale = FALSE)
    meanX = apply(Xtrain, 2, mean)
    meanY = apply(Ytrain, 2, mean)
    n = dim(X)[1]
    p = dim(X)[2]
    m = dim(Y)[2]
    if (is.null(ncomp)){
       ncomp = min(n,p)
    }
    if (weight == TRUE){
       q = matrix(, n, 1)
       for(i in 1:n){
           q[i,] = beta * (1 - beta)^(i - 1)
       }
       weig = diag(rev(q))
       X = weig %*% X
       Y = weig %*% as.matrix(Y)
    }
    S = crossprod(X, Y)
    RR = matrix(0, ncol = max(ncomp), nrow = p)
    PP = matrix(0, ncol = max(ncomp), nrow = p)
    QQ = matrix(0, ncol = max(ncomp), nrow = m)
    TT = matrix(0, ncol = max(ncomp), nrow = n)
    VV = matrix(0, ncol = max(ncomp), nrow = p)
    UU = matrix(0, ncol = max(ncomp), nrow = n)
    B = array(0, c(dim(X)[2], dim(Y)[2], length(ncomp)))
    if (!is.null(Xtest)) {
        if (is.vector(Xtest)) {
            Xtest = matrix(Xtest, 1, length(Xtest))
        }
        Ypred = array(0, c(dim(Xtest)[1], m, length(ncomp)))
    }
    for (a in 1:max(ncomp)) {
        qq = svd(S)$v[, 1]
        rr = S %*% qq
        tt = scale(X %*% rr, scale = FALSE)
        tnorm = sqrt(sum(tt * tt))
        tt = tt / tnorm
        rr = rr / tnorm
        pp = crossprod(X, tt)
        qq = crossprod(Y, tt)
        uu = Y %*% qq
        vv = pp
        if (a > 1) {
            vv = vv - VV %*% crossprod(VV, pp)
            uu = uu - TT %*% crossprod(TT, uu)
        }
        vv = vv / sqrt(sum(vv * vv))
        S = S - vv %*% crossprod(vv, S)
        RR[, a] = rr
        TT[, a] = tt
        PP[, a] = pp
        QQ[, a] = qq
        VV[, a] = vv
        UU[, a] = uu
        if (!is.na(i <- match(a, ncomp))) {
            B[, , i] = RR[, 1:a, drop = FALSE] %*% t(QQ[, 1:a, 
                drop = FALSE])
            if (!is.null(Xtest)) {
                Xtest = scale(Xtest, scale = FALSE, center = meanX)
                Ypred[, , i] = Xtest %*% B[, , i] + meanY
            }
        }
    }
    if (length(ncomp) == 1) {
        B = B[, , 1]
    }
    if (!is.null(Xtest)) 
        if (weight==TRUE){
            list(B = B, Ypred = Ypred, P = PP, Q = QQ, T = TT, R = RR, weight = q)
        }
        else { 
            list(B = B, Ypred = Ypred, P = PP, Q = QQ, T = TT, R = RR)
        } 
    else 
        if (weight==TRUE){
            list(B = B, P = PP, Q = QQ, T = TT, R = RR, weight = q)
        }
        else {
            list(B = B, P = PP, Q = QQ, T = TT, R = RR)
        }
}