Ngme2 random walk model

In this vignette, we will briefly introduce the random walk model in ngme2.

Description

In ngme2, we currently support random walk model of the first-order and the second order.

Remeber in ngme2, all the models has the same form: $$K \mathbf{W} = \boldsymbol{\epsilon},$$ where $K$ is the operator matrix, $\epsilon$ is the random noise (maybe non-Gaussian), $\mathbf{W}$ is the random vector we want to model.

The first-order random walk (RW(1))

The first-order random walk is constructed assuming independent increments given data $\mathbf{x}$ (use as location, sorted) of length $n$:

1. When $\epsilon$ follows the NIG distribution: $$\begin{align} \Delta w_i &= -\mu + \mu V_i + \sigma \sqrt{V_i} Z , \; n=1, \dots, n-1 \\ V_i &\sim IG(\nu, \nu (\Delta x_i)^2), \\ Z_i &\sim N(0, 1), \end{align}$$

2. When $\epsilon$ follows the normal distribution: $$\begin{align} \Delta w_i &= \epsilon_i , \; n=1, \dots, n-1 \\ \epsilon_i &\sim N(0, \sigma^2 \Delta x_i), \\ \end{align}$$

where $\Delta w_i := w_{i+1} - w_{i}$, $\Delta x_i := x_{i+1} - x_{i}$.

The operator matrix $K$ of dimension $(n-1 \times n)$ is

$$K = \begin{bmatrix} -1 & 1 \\ & -1 & 1 \\ & & \ddots & \ddots \\ & & & -1 & 1 \end{bmatrix}.$$

We also provide the special case of circular random walk, which the 1st element and n-th element is connected. In the circular RW(1) case, the operator matrix $K$ of dimension $(n-1 \times n-1)$ is

$$K = \begin{bmatrix} -1 & 1 \\ & \ddots & \ddots \\ & & -1 & 1 \\ 1 & & & -1 \end{bmatrix}.$$

The second-order random walk (RW(2))

Similarily, the second-order random walk is constructed assuming the second order difference is independent:

1. When $\epsilon$ follows the NIG distribution: $$\begin{align} \Delta^2 w_i &= -\mu + \mu V_i + \sigma \sqrt{V_i} Z , \; n=1, \dots, n-2 \\ V_i &\sim IG(\nu, \nu (\Delta x_i)^2), \\ Z_i &\sim N(0, 1), \end{align}$$

2. When $\epsilon$ follows the normal distribution: $$\begin{align} \Delta^2 w_i &= \epsilon_i , \; n=1, \dots, n-2 \\ \epsilon_i &\sim N(0, \sigma^2 \Delta x_i), \\ \end{align}$$

where $\Delta^2 w_i := w_{i+2} - 2w_{i+1} - w_{i}$, $\Delta x_i := x_{i+1} - x_{i}$.

The operator matrix $K$ of dimension $(n-2 \times n)$ is

$$K = \begin{bmatrix} 1 & -2 & 1 \\ & 1 & -2 & 1 \\ & & \ddots & \ddots & \ddots \\ & & & 1 & -2 & 1 \end{bmatrix}.$$

In the circular RW(2) case, the operator matrix $K$ of dimension $(n-2 \times n-2)$ is

$$K = \begin{bmatrix} 1 & -2 & 1 \\ & 1 & -2 & 1 \\ & & \ddots & \ddots & \ddots \\ 1 & & & 1 & -2 \\ -2 & 1 & & & 1 \end{bmatrix}.$$

Usage

Use the f(model = "rw1") or f(model = "rw2") (in formula) to specify a random walk model.

library(ngme2)
#> This is ngme2 of version 0.7.0
#> - See our homepage: https://davidbolin.github.io/ngme2 for more details.
#> 
#> Attaching package: 'ngme2'
#> The following object is masked from 'package:stats':
#> 
#>     ar
set.seed(16)
m1 <- f(rexp(5), model="rw1", noise = noise_normal())
m1$operator$K
#> 5 x 5 sparse Matrix of class "dgCMatrix"
#>                                                            
#> [1,]  0.02295832  0.06577409  0.1130829  0.730911 0.6606439
#> [2,] -1.00000000  1.00000000  .          .        .        
#> [3,]  .          -1.00000000  1.0000000  .        .        
#> [4,]  .           .          -1.0000000  1.000000 .        
#> [5,]  .           .           .         -1.000000 1.0000000

m2 <- f(rnorm(6), model="rw2", cyclic = TRUE)
m2$operator$K
#> 7 x 7 sparse Matrix of class "dgCMatrix"
#>                                                                        
#> [1,] .           0.02159833  0.6093765  0.6305594  0.0630014  0.3724705
#> [2,] 0.02159833  1.00000000 -2.0000000  1.0000000  .          .        
#> [3,] 0.60937647  .           1.0000000 -2.0000000  1.0000000  .        
#> [4,] 0.63055937  .           .          1.0000000 -2.0000000  1.0000000
#> [5,] 0.06300140  .           .          .          1.0000000 -2.0000000
#> [6,] 0.37247049  1.00000000  .          .          .          1.0000000
#> [7,] 0.35225032 -2.00000000  1.0000000  .          .          .        
#>                
#> [1,]  0.3522503
#> [2,]  .        
#> [3,]  .        
#> [4,]  .        
#> [5,]  1.0000000
#> [6,] -2.0000000
#> [7,]  1.0000000

Simulation

Doing the simulation is simple, just pass the corresponding model into simulate function.

simulate(m1)
#>         sim_1
#> 1  0.42732713
#> 2  0.13792977
#> 3 -0.01727820
#> 4 -0.05194693
#> 5  0.24731142