knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library(ggplot2) library(SMFilter)
The package implements the filtering algorithms for the state-space models on the Stiefel manifold. It also implements sampling algorithms for uniform, vector Langevin-Bingham and matrix Langevin-Bingham distributions on the Stiefel manifold. You can also find the package on CRAN, see
and the corresponding paper
State-Space Models on the Stiefel Manifold with a New Approach to Nonlinear Filtering
You can either install the stable version from CRAN
install.packages("SMFilter")
or install the development version from GitHub
devtools::install_github("yukai-yang/SMFilter")
provided that the package "devtools" has been installed beforehand.
After installing the package, you need to load (attach better say) it by running the code
library(SMFilter)
You can first check the information and the current version number by running
version()
Then you can take a look at all the available functions and data in the package
ls( grep("SMFilter", search()) )
For details, see
?SimModel1
First we can use the package to sample from the type one model. To this end, we shall initialize by running
set.seed(1) # control the seed iT = 100 # sample size ip = 2 # dimension of the dependent variable ir = 1 # rank number iqx = 3 # dimension of the independent variable x_t iqz=0 # dimension of the independent variable z_t ik = 0 # lag length method='max_3' # the optimization methond to use, for details, see FilterModel1 Omega = diag(ip)*.1 # covariance of the errors vD = 50 # diagonal of the D matrix
Then we initialize the data and some other parameters
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) if(ip*ik+iqz==0) mB=NULL else mB = matrix(c(runif_sm(num=1,ip=(ip*ik+iqz)*ip,ir=1)), ip, ip*ik+iqz)
Then we can simulate from the model
ret = SimModel1(iT=iT, mX=mX, mZ=mZ, mY=mY, alpha_0=alpha_0, beta=beta, mB=mB, vD=vD, Omega=Omega)
Have a look at the simulated data
matplot(ret$dData[,1:ip], type="l", ylab="simulated data")
Then let's apply the filtering algorithm on the data
fil = FilterModel1(mY=as.matrix(ret$dData[,1:ip]), mX=mX, mZ=mZ, beta=beta, mB=mB, Omega=Omega, vD=vD, U0=alpha_0, method=method)
Then we compare the filtered modal orientations with the true ones in terms of the Frobenius distance.
fil = fil[2:(iT+1),,,drop=F] # remove the initial value ra = ret$aAlpha # get the true ones # define a function to compute the distances ftmp <- function(ix){ mx1 = matrix(fil[ix,,],ip,ir) mx2 = matrix(ra[ix,,],ip,ir) return(FDist2(mx1,mx2)) } # plot the distances ggplot() + geom_point(aes(x=1:iT,y=sapply(1:iT,ftmp)/4/ir)) + ylim(0, 1) + labs(x=paste0("t= 1, ...,",iT), y="normalized d( \u03b1, U )")
For details, see
?SimModel2
Again, we start with sampling. We initialize the parameters
iT = 100 ip = 2 ir = 1 iqx = 4 iqz=0 ik = 0 Omega = diag(ip)*.1 vD = 50
Then we initialize the data and some other parameters
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) if(ip*ik+iqz==0) mB=NULL else mB = matrix(c(runif_sm(num=1,ip=(ip*ik+iqz)*ip,ir=1)), ip, ip*ik+iqz)
Then we can simulate from the model
ret = SimModel2(iT=iT, mX=mX, mZ=mZ, mY=mY, alpha=alpha, beta_0=beta_0, mB=mB, vD=vD)
And then have a look at the simulated data
matplot(ret$dData[,1:ip], type="l",ylab="simulated data")
Apply the filtering algorithm on the data
fil = FilterModel2(mY=as.matrix(ret$dData[,1:ip]), mX=mX, mZ=mZ, alpha=alpha, mB=mB, Omega=Omega, vD=vD, U0=beta_0, method=method)
Then we compare the filtered modal orientations with the true ones in terms of the Frobenius distance.
fil = fil[2:(iT+1),,,drop=F] # remove the initial value ra = ret$aBeta ftmp <- function(ix){ mx1 = matrix(fil[ix,,],iqx,ir) mx2 = matrix(ra[ix,,],iqx,ir) return(FDist2(mx1,mx2)) } # plot the distances ggplot() + geom_point(aes(x=1:iT,y=sapply(1:iT,ftmp)/4/ir)) + ylim(0, 1) + labs(x=paste0("t= 1, ...,",iT), y="normalized d( \u03b2, U )")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.