# R/filtering.R In SMFilter: Filtering Algorithms for the State Space Models on the Stiefel Manifold

#### Documented in FilterModel1FilterModel2

#################################################################################
## filtering algorithms
## Nov 2018
#################################################################################

#' Filtering algorithm for the type one model.
#'
#' This function implements the filtering algorithm for the type one model.
#' See Details part below.
#'
#' The type one model on Stiefel manifold takes the form:
#' \deqn{\boldsymbol{y}_t \quad = \quad \boldsymbol{\alpha}_t \boldsymbol{\beta} ' \boldsymbol{x}_t + \boldsymbol{B} \boldsymbol{z}_t + \boldsymbol{\varepsilon}_t}
#' \deqn{\boldsymbol{\alpha}_{t+1} | \boldsymbol{\alpha}_{t} \quad \sim \quad ML (p, r, \boldsymbol{\alpha}_{t} \boldsymbol{D})}
#' where \eqn{\boldsymbol{y}_t} is a \eqn{p}-vector of the dependent variable,
#' \eqn{\boldsymbol{x}_t} and \eqn{\boldsymbol{z}_t} are explanatory variables wit dimension \eqn{q_1} and \eqn{q_2},
#' \eqn{\boldsymbol{x}_t} and \eqn{\boldsymbol{z}_t} have no overlap,
#' matrix \eqn{\boldsymbol{B}} is the coefficients for \eqn{\boldsymbol{z}_t},
#' \eqn{\boldsymbol{\varepsilon}_t} is the error vector.
#'
#' The matrices \eqn{\boldsymbol{\alpha}_t} and \eqn{\boldsymbol{\beta}} have dimensions \eqn{p \times r} and \eqn{q_1 \times r}, respectively.
#' Note that \eqn{r} is strictly smaller than both \eqn{p} and \eqn{q_1}.
#' \eqn{\boldsymbol{\alpha}_t} and \eqn{\boldsymbol{\beta}} are both non-singular matrices.
#' \eqn{\boldsymbol{\alpha}_t} is time-varying while \eqn{\boldsymbol{\beta}} is time-invariant.
#'
#' Furthermore, \eqn{\boldsymbol{\alpha}_t} fulfills the condition \eqn{\boldsymbol{\alpha}_t' \boldsymbol{\alpha}_t = \boldsymbol{I}_r},
#' and therefor it evolves on the Stiefel manifold.
#'
#' \eqn{ML (p, r, \boldsymbol{\alpha}_{t} \boldsymbol{D})} denotes the Matrix Langevin distribution or matrix von Mises-Fisher distribution on the Stiefel manifold.
#' Its density function takes the form
#' \deqn{f(\boldsymbol{\alpha_{t+1}} ) = \frac{ \mathrm{etr} \left\{ \boldsymbol{D} \boldsymbol{\alpha}_{t}' \boldsymbol{\alpha_{t+1}} \right\} }{ _{0}F_1 (\frac{p}{2}; \frac{1}{4}\boldsymbol{D}^2 ) }}
#' where \eqn{\mathrm{etr}} denotes \eqn{\mathrm{exp}(\mathrm{tr}())},
#' and \eqn{_{0}F_1 (\frac{p}{2}; \frac{1}{4}\boldsymbol{D}^2 )} is the (0,1)-type hypergeometric function for matrix.
#'
#'
#' @param mY the matrix containing Y_t with dimension \eqn{T \times p}.
#' @param mX the matrix containing X_t with dimension \eqn{T \times q_1}.
#' @param mZ the matrix containing Z_t with dimension \eqn{T \times q_2}.
#' @param beta the \eqn{\beta} matrix.
#' @param mB the coefficient matrix \eqn{\boldsymbol{B}} before \code{mZ} with dimension \eqn{p \times q_2}.
#' @param Omega covariance matrix of the errors.
#' @param vD vector of the diagonals of \eqn{D}.
#' @param U0 initial value of the alpha sequence.
#' @param method a string representing the optimization method from c('max_1','max_2','max_3','min_1','min_2').
#'
#' @return an array \code{aAlpha} containing the modal orientations of alpha in the prediction step.
#'
#' @author Yukai Yang, \email{yukai.yang@@statistik.uu.se}
#' @keywords filtering
#'
#' @examples
#'
#' iT = 50
#' ip = 2
#' ir = 1
#' iqx = 4
#' iqz=0
#' ik = 0
#' Omega = diag(ip)*.1
#'
#' if(iqx==0) mX=NULL else mX = matrix(rnorm(iT*iqx),iT, iqx)
#' if(iqz==0) mZ=NULL else mZ = matrix(rnorm(iT*iqz),iT, iqz)
#' if(ik==0) mY=NULL else mY = matrix(0, ik, ip)
#'
#' alpha_0 = matrix(c(runif_sm(num=1,ip=ip,ir=ir)), ip, ir)
#' beta = matrix(c(runif_sm(num=1,ip=ip*ik+iqx,ir=ir)), ip*ik+iqx, ir)
#' mB=NULL
#' vD = 100
#'
#' ret = SimModel1(iT=iT, mX=mX, mZ=mZ, mY=mY, alpha_0=alpha_0, beta=beta, mB=mB, vD=vD, Omega=Omega)
#' mYY=as.matrix(ret$dData[,1:ip]) #' fil = FilterModel1(mY=mYY, mX=mX, mZ=mZ, beta=beta, mB=mB, Omega=Omega, vD=vD, U0=alpha_0) #' #' @export FilterModel1 <- function(mY, mX, mZ, beta, mB=NULL, Omega, vD, U0, method='max_1') { iT = nrow(mY) ip = ncol(mY) ir = ncol(beta) aAlpha = array(0,dim=c(iT+1,ip,ir)) mJ = chol2inv(chol(Omega)) mD = diag(vD,ir,ir) if(is.null(mB)){ mBZ = matrix(0,iT,ip) }else{ mBZ = NULL for(iter in 1:iT) mBZ = rbind(mBZ, mB%*%mZ[iter,]) } if(length(method) > 1) stop(simpleError("Don't choose multiple methods!")) opt = switch(match(method,c('max_1','max_2','max_3','min_1','min_2')), MaxFx1, MaxFx2, MaxFx3, MinLine, MinLineOpt) if(is.null(opt)) stop(simpleError("The method does not match from c('max_1','max_2','max_3','min_1','min_2')!")) mU = U0 aAlpha[1,,] = mU for(iter in 1:iT){ mH = -.5 * t(beta)%*%c(mX[iter,])%*%t(mX[iter,])%*%beta mC = mU%*%mD + mJ%*% c(mY[iter,]-mBZ[iter,])%*%t(mX[iter,])%*%beta tmp = svd(mC) tmp = tmp$u%*%t(tmp$v) # maximize and update mU mU = opt(mX=tmp,mA=mH,mB=mJ,mC=mC)$mX

aAlpha[iter+1,,] = mU
}
return(aAlpha)
}

#' Filtering algorithm for the type two model.
#'
#' This function implements the filtering algorithm for the type two model.
#' See Details part below.
#'
#' The type two model on Stiefel manifold takes the form:
#' \deqn{\boldsymbol{y}_t \quad = \quad \boldsymbol{\alpha} \boldsymbol{\beta}_t ' \boldsymbol{x}_t + \boldsymbol{B}' \boldsymbol{z}_t + \boldsymbol{\varepsilon}_t}
#' \deqn{\boldsymbol{\beta}_{t+1} | \boldsymbol{\beta}_{t} \quad \sim \quad ML (q_1, r, \boldsymbol{\beta}_{t} \boldsymbol{D})}
#' where \eqn{\boldsymbol{y}_t} is a \eqn{p}-vector of the dependent variable,
#' \eqn{\boldsymbol{x}_t} and \eqn{\boldsymbol{z}_t} are explanatory variables wit dimension \eqn{q_1} and \eqn{q_2},
#' \eqn{\boldsymbol{x}_t} and \eqn{\boldsymbol{z}_t} have no overlap,
#' matrix \eqn{\boldsymbol{B}} is the coefficients for \eqn{\boldsymbol{z}_t},
#' \eqn{\boldsymbol{\varepsilon}_t} is the error vector.
#'
#' The matrices \eqn{\boldsymbol{\alpha}} and \eqn{\boldsymbol{\beta}_t} have dimensions \eqn{p \times r} and \eqn{q_1 \times r}, respectively.
#' Note that \eqn{r} is strictly smaller than both \eqn{p} and \eqn{q_1}.
#' \eqn{\boldsymbol{\alpha}} and \eqn{\boldsymbol{\beta}_t} are both non-singular matrices.
#' \eqn{\boldsymbol{\beta}_t} is time-varying while \eqn{\boldsymbol{\alpha}} is time-invariant.
#'
#' Furthermore, \eqn{\boldsymbol{\beta}_t} fulfills the condition \eqn{\boldsymbol{\beta}_t' \boldsymbol{\beta}_t = \boldsymbol{I}_r},
#' and therefor it evolves on the Stiefel manifold.
#'
#' \eqn{ML (p, r, \boldsymbol{\beta}_t \boldsymbol{D})} denotes the Matrix Langevin distribution or matrix von Mises-Fisher distribution on the Stiefel manifold.
#' Its density function takes the form
#' \deqn{f(\boldsymbol{\beta_{t+1}} ) = \frac{ \mathrm{etr} \left\{ \boldsymbol{D} \boldsymbol{\beta}_{t}' \boldsymbol{\beta_{t+1}} \right\} }{ _{0}F_1 (\frac{p}{2}; \frac{1}{4}\boldsymbol{D}^2 ) }}
#' where \eqn{\mathrm{etr}} denotes \eqn{\mathrm{exp}(\mathrm{tr}())},
#' and \eqn{_{0}F_1 (\frac{p}{2}; \frac{1}{4}\boldsymbol{D}^2 )} is the (0,1)-type hypergeometric function for matrix.
#'
#'
#' @param mY the matrix containing Y_t with dimension \eqn{T \times p}.
#' @param mX the matrix containing X_t with dimension \eqn{T \times q_1}.
#' @param mZ the matrix containing Z_t with dimension \eqn{T \times q_2}.
#' @param alpha the \eqn{\alpha} matrix.
#' @param mB the coefficient matrix \eqn{\boldsymbol{B}} before \code{mZ} with dimension \eqn{p \times q_2}.
#' @param Omega covariance matrix of the errors.
#' @param vD vector of the diagonals of \eqn{D}.
#' @param U0 initial value of the alpha sequence.
#' @param method a string representing the optimization method from c('max_1','max_2','max_3','min_1','min_2').
#'
#' @return an array \code{aAlpha} containing the modal orientations of alpha in the prediction step.
#'
#' @author Yukai Yang, \email{yukai.yang@@statistik.uu.se}
#' @keywords filtering
#'
#' @examples
#'
#' iT = 50
#' ip = 2
#' ir = 1
#' iqx = 4
#' iqz=0
#' ik = 0
#' Omega = diag(ip)*.1
#'
#' if(iqx==0) mX=NULL else mX = matrix(rnorm(iT*iqx),iT, iqx)
#' if(iqz==0) mZ=NULL else mZ = matrix(rnorm(iT*iqz),iT, iqz)
#' if(ik==0) mY=NULL else mY = matrix(0, ik, ip)
#'
#' alpha = matrix(c(runif_sm(num=1,ip=ip,ir=ir)), ip, ir)
#' beta_0 = matrix(c(runif_sm(num=1,ip=ip*ik+iqx,ir=ir)), ip*ik+iqx, ir)
#' mB=NULL
#' vD = 100
#'
#' ret = SimModel2(iT=iT, mX=mX, mZ=mZ, mY=mY, alpha=alpha, beta_0=beta_0, mB=mB, vD=vD)
#' mYY=as.matrix(ret$dData[,1:ip]) #' fil = FilterModel2(mY=mYY, mX=mX, mZ=mZ, alpha=alpha, mB=mB, Omega=Omega, vD=vD, U0=beta_0) #' #' @export FilterModel2 <- function(mY, mX, mZ, alpha, mB=NULL, Omega, vD, U0, method='max_1') { iT = nrow(mY) ip = ncol(mY) iqxx = ncol(mX) ir = ncol(alpha) aBeta = array(0,dim=c(iT+1,iqxx,ir)) iOmg = chol2inv(chol(Omega)) mH = -.5*t(alpha)%*%iOmg%*%alpha mU = U0 aBeta[1,,] = mU mD = diag(vD,ir,ir) if(is.null(mB)){ mBZ = matrix(0,iT,ip) }else{ mBZ = NULL for(iter in 1:iT) mBZ = rbind(mBZ, mB%*%mZ[iter,]) } if(length(method) > 1) stop(simpleError("Don't choose multiple methods!")) opt = switch(match(method,c('max_1','max_2','max_3','min_1','min_2')), MaxFx1, MaxFx2, MaxFx3, MinLine, MinLineOpt) if(is.null(opt)) stop(simpleError("The method does not match from c('max_1','max_2','max_3','min_1','min_2')!")) for(iter in 1:iT){ mJ = mX[iter,] %*% t(mX[iter,]) mC = mU%*%mD + mX[iter,]%*%t(mY[iter,]-mBZ[iter,])%*%iOmg%*%alpha # maximize and update mU mU = opt(mX=mU,mA=mH,mB=mJ,mC=mC)$mX

aBeta[iter+1,,] = mU
}
return(aBeta)
}


## Try the SMFilter package in your browser

Any scripts or data that you put into this service are public.

SMFilter documentation built on May 1, 2019, 8:01 p.m.