---
title: "SMFilter: Filtering Algorithms for the State-Space models on the Stiefel Manifold"
author: "Yukai Yang"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{README}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(ggplot2)
library(SMFilter)
```

# SMFilter version 1.0.3 (Red Filter)


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

[SMFilter@CRAN](https://CRAN.R-project.org/package=SMFilter)

and the corresponding paper

[State-Space Models on the Stiefel Manifold with a New Approach to Nonlinear Filtering](https://www.mdpi.com/2225-1146/6/4/48) 


## How to install

You can either install the stable version from CRAN
```{r install1, eval=F}
install.packages("SMFilter")
```
or install the development version from GitHub
```{r install2, eval=F}
devtools::install_github("yukai-yang/SMFilter")
```
provided that the package "devtools" has been installed beforehand.

## Example

After installing the package, you need to load (attach better say) it by running the code
```{r attach}
library(SMFilter)
```

You can first check the information and the current version number by running
```{r version}
version()
```

Then you can take a look at all the available functions and data in the package
```{r contents}
ls( grep("SMFilter", search()) ) 
```


### Type one model

For details, see

```{r type1, eval=F}
?SimModel1
```


First we can use the package to sample from the type one model. To this end, we shall initialize by running

```{r init11}
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

```{r init12}
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

```{r sim1}
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

```{r dat1}
matplot(ret$dData[,1:ip], type="l", ylab="simulated data")
```

Then let's apply the filtering algorithm on the data

```{r filter1}
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.

```{r dist1, echo=F}
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 )")
```



### Type two model

For details, see

```{r type2, eval=F}
?SimModel2
```

Again, we start with sampling. We initialize the parameters

```{r init21}
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

```{r init22}
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

```{r sim2}
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

```{r dat2}
matplot(ret$dData[,1:ip], type="l",ylab="simulated data")
```

Apply the filtering algorithm on the data

```{r filter2}
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.

```{r dist, echo=F}
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 )")
```