In jlin-vt/SML: Statistical Machine Learning

knitr::opts_chunk$set(echo = TRUE)
library(SML)
library(mvtnorm)
library(pscl)
library(latex2exp)
library(ggplot2)

1. Example: Sampling for univariate random variate using its CDF

Goal: sample from $f(x) = (x - 1) ^ 2 * exp{-(x ^ 3 / 3 - 2 * x ^ 2 / 2 + x)}$.
$f(x)$: density function.
$F(x)$: cdf of $f(x)$.
$F^{-1}(y)$ = inverse cdf of $f(x)$.

Define a target pdf.

f <- function(x) {((x - 1) ** 2) * exp(-(x ** 3 / 3 - 2 * x ** 2 / 2 + x))}

Load high-order function from the package.

F = CDF(f)
F.inv = InverseCDF(F(x), l=0, u=10)

Generate regular sequences.

x <- seq(0, 5, length.out = 1000)
y <- seq(0, 1, length.out = 1000)

par(mfrow=c(1, 3))
plot(x, f(x), type = "l", main = "f(x)")
plot(x, F(x), type = "l", main = "CDF of f(x)")
plot(y, F.inv(y), type = "l", main = "Inverse CDF of f(x)")

Generate a random sample and make a plot.

X <- runif(1000, 0, 1)  
Z <- F.inv(X)

par(mfrow = c(1, 2))
plot(x, f(x), type = "l",main = "Density function")
hist(Z, breaks = 20, xlim = c(0, 5))

2. Example: Rejection sampling from target density

Target density: $f(\theta) = (2 / pi)^{\frac{1}{2}} * \exp(-\theta^ 2 / 2)$.
Proposal distribution: $g(\theta) = \exp(-\theta)$

Define target and proposal

theta.list <- seq(0, 10, length.out = 100)
f <- function(x) sqrt(2 / pi) * exp(-x ** 2/ 2)
g <- function(x) exp(-x)
M <- 1.4

Plot $f$ and $g$:

sample = RejectionSampler(theta.list, f, g, M)
plot(theta.list, f(theta.list), type = "l", lwd = 2)
lines(theta.list, g(theta.list), col = "blue", lwd = 3)

hist(sample, breaks = 20, probability = TRUE, add = TRUE, col = "orange")

3. Example: SMC method using SIR algorithm in linear Gaussian state-space model

$x$: latent state.
$y$: observed state.
$x_1 \sim N(0, 1)$.
$x_n = \phi x_{n-1} + v_n, v_n \sim N(0,s_v^2)$.
$y_n = x_n + w_n, w_n \sim N(0,s_w^2)$.

Generate data.

set.seed(100)
T <- 100
S.v <- 1
S.w <- 1 
phi <- 0.95
x <- rep(NA, T)   
y <- rep(NA, T)   
x[1] <- rnorm(1, mean = 0, sd = 1)

for (i in 2:T){
  x[i] <- phi*x[i-1] + rnorm(1, mean = 0, sd = S.v)
}
y <- x + rnorm(T, mean = 0, sd = S.w)

Run sampling.

res = SirLinearGauss(y, T, N, S.v, S.w, phi)
mu.Xs = res[["mu.Xs"]]  
var.Xs = res[["var.Xs"]]

Plot.

plot(1:T, x, type = 'b', pch = 3, col = 'grey', main = 'SIR Solution', ylim = c(-4, 5))
points(1:T, y, pch = 16, col = 'blue')
lines(1:T, y,  col = 'blue')
lines(mu.Xs - 2 * sqrt(var.Xs), lty = 2, col = 'black')
lines(mu.Xs + 2 * sqrt(var.Xs), lty = 2, col = 'black')
lines(mu.Xs, col = "red", lwd = 2)
points(mu.Xs, col = "red", lwd = 2, pch = 16)
legend(10, -2.4, c('latent', 'observed', 'post.mean', 'post.CI'),
       pch = c(3, 16, 16, NA), col = c('grey', 'blue', 'red', 'black'), lty = c(1, 1, 1, 2))

4. Example: SMC method using SIR algorithm in Stochastic volatility model.

$x$: latent state.
$y$: observed state.
$x_1 \sim N(0, \frac{\sigma^2}{1-\sigma^2})$
$x_n = \phi x_{n-1} + v_n, v_n \sim N(0, \sigma^2), |sigma^2|<1$
$y_n = exp{\gamma + x_n}w_n, w_n \sim N(0, 1)$

Initial value:

Time <- 100
sigma <- 1
gamma <- 1 
phi <- 0.9
x <- rep(NA, Time)  
y <- rep(NA, Time)  
x[1] <- rnorm(1, mean = 0, sd = sigma / sqrt(1 - phi^2))
for (i in 2:Time){
  x[i] <- phi * x[i-1] + rnorm(1, mean = 0, sd = sigma)
}
y <- exp(gamma + x) * rnorm(Time, mean = 0, sd = 1)

Plot the latent sequence and the observed sequence:

plot(1:Time, x, type = 'b', pch = 3, col = 'green', xlab = "Time", ylab = NULL)
points(1:Time, y, pch = 16, col = 'blue')
lines(1:Time, y,  col = 'blue')
legend(-1, -3, c('latent', 'observed'), col = c('green','blue'), pch = c(3, 16), lty = c(1, 1))

Run sampling.

res = SirStochasticVolatility(y, T, N, gamma, sigma, phi)
mu.Xs = res[["mu.Xs"]]  
var.Xs = res[["var.Xs"]]

Plot:

plot(1:Time, x, type = 'b', pch = 3, col = 'grey', main = 'SIR Solution', ylim = c(-10, 20), xlab = "Time")
points(1:Time, y, pch = 16, col = 'blue')
lines(1:Time, y,  col = 'blue')
lines(mu.Xs - 2 * sqrt(var.Xs), lty = 2, col = 'black') 
lines(mu.Xs + 2 * sqrt(var.Xs), lty = 2, col = 'black')
lines(mu.Xs, col = "red",lwd = 2)
points(mu.Xs, col = "red",lwd = 2, pch = 16)

legend(0, -5, c('latent', 'observed', 'post.mean', 'post.CI'),
       pch = c(3, 16, 16, NA), col = c('grey', 'blue', 'red', 'black'), lty = c(1, 1, 1, 2))