Expectation-Maximization

Description

expectation_maximization is a method to choose hyper-parameters of the linear Gaussian State-Space Model with time-invariant variances relying on the Expectation-Maximization algorithm.

Usage

expectation_maximization(
  X,
  y,
  n_iter,
  Q_init,
  sig_init = 1,
  verbose = 1000,
  lambda = 10^-9,
  mode_diag = FALSE,
  p1 = 0
)

Arguments

Details

E-step is realized through recursive Kalman formulae (filtering then smoothing).
M-step is the maximization of the expected complete likelihood with respect to the hyper-parameters.
We only have the guarantee of convergence to a LOCAL optimum. We fix P1 = p1 I (by default p1 = 0). We optimize theta1, sig, Q.

Value

a list containing values for P,theta,sig,Q, and two vectors DIFF, LOGLIK assessing the convergence of the algorithm.

Examples

set.seed(1)
### Simulate data
n <- 100
d <- 5
Q <- diag(c(0,0,0.25,0.25,0.25))
sig <- 1

X <- cbind(matrix(rnorm((d-1)*n,sd=1),n,d-1),1)
theta <- matrix(rnorm(d), d, 1)
theta_arr <- matrix(0, n, d)
for (t in 1:n) {
  theta_arr[t,] <- theta
  theta <- theta + matrix(mvtnorm::rmvnorm(1,matrix(0,d,1),Q),d,1)
}
y <- rowSums(X * theta_arr) + rnorm(n, sd=sig)

l <- viking::expectation_maximization(X, y, 50, diag(d), verbose=10)
print(l$Q)
print(l$sig)

Iterative Grid Search

Description

iterative_grid_search is an iterative method to choose hyper-parameters of the linear Gaussian State-Space Model with time-invariant variances.

Usage

iterative_grid_search(
  X,
  y,
  q_list,
  Q_init = NULL,
  max_iter = 0,
  delay = 1,
  use = NULL,
  restrict = NULL,
  mode = "gaussian",
  p1 = 0,
  ncores = 1,
  train_theta1 = NULL,
  train_Q = NULL,
  verbose = TRUE
)

Arguments

Details

We restrict ourselves to a diagonal matrix Q and we optimize Q / sig^2 on a grid. Each diagonal coefficient is assumed to belong to a pre-defined q_list.
We maximize the log-likelihood on that region of search in an iterative fashion. At each step, we change the diagonal coefficient improving the most the log-likelihood. We stop when there is no possible improvement. This doesn't guarantee an optimal point on the grid, but the computational time is much lower.

Value

a list containing values for P,theta,sig,Q, as well as LOGLIK, the evolution of the log-likelihood during the search.

Examples

set.seed(1)
### Simulate data
n <- 100
d <- 5
Q <- diag(c(0,0,0.25,0.25,0.25))
sig <- 1

X <- cbind(matrix(rnorm((d-1)*n,sd=1),n,d-1),1)
theta <- matrix(rnorm(d), d, 1)
theta_arr <- matrix(0, n, d)
for (t in 1:n) {
  theta_arr[t,] <- theta
  theta <- theta + matrix(mvtnorm::rmvnorm(1,matrix(0,d,1),Q),d,1)
}
y <- rowSums(X * theta_arr) + rnorm(n, sd=sig)

l <- viking::iterative_grid_search(X, y, seq(0,1,0.25))
print(l$Q)
print(l$sig)

Kalman Filtering

Description

Compute the filtered estimation of the parameters theta and P.

Usage

kalman_filtering(X, y, theta1, P1, Q = 0, sig = 1)

Arguments

Value

a list containing theta_arr and P_arr, the filtered estimation of the parameters theta and P.

Kalman Smoothing

Description

Compute the smoothed estimation of the parameters theta and P.

Usage

kalman_smoothing(X, y, theta1, P1, Q = 0, sig = 1)

Arguments

Value

a list containing theta_arr and P_arr, the smoothed estimation of the parameters theta and P.

Log-likelihood

Description

loglik computes the log-likelihood of a state-space model of specified Q/sig^2, P1/sig^2, theta1.

Usage

loglik(X, y, Qstar, use, p1, train_theta1, train_Q, mode = "gaussian")

Arguments

Value

a numeric value for the log-likelihood.

Plot a statespace object

Description

plot.statespace displays different graphs expressing the behavior of the state-space model:
1. Evolution of the Bias: rolling version of the error of the model.
2. Evolution of the RMSE: root-mean-square-error computed on a rolling window.
3. State Evolution: time-varying state coefficients, subtracted of the initial state vector.
4. Normal Q-Q Plot: we check if the observation follows the Gaussian distribution of estimated mean and variance. To that end, we display a Q-Q plot of the residual divided by the estimated standard deviation, against the standard normal distribution.

Usage

## S3 method for class 'statespace'
plot(x, pause = FALSE, window_size = 7, date = NULL, sel = NULL, ...)

Arguments

Value

No return value, called to display plots.

Predict using a statespace object

Description

predict.statespace makes a prediction for a statespace object, in the offline or online setting.

Usage

## S3 method for class 'statespace'
predict(
  object,
  newX,
  newy = NULL,
  online = TRUE,
  compute_smooth = FALSE,
  type = c("mean", "proba", "model"),
  ...
)

Arguments

Value

Depending on the type specified, the result is
- a vector of mean forecast if type='mean' - a list of two vectors, mean forecast and standard deviations if type='proba' - a statespace object if type='model'

Select time-invariant variances of a State-Space Model

Description

select_Kalman_variances is a function to choose hyper-parameters of the linear Gaussian State-Space Model with time-invariant variances. It relies on the functions iterative_grid_search and expectation_maximization.

Usage

select_Kalman_variances(ssm, X, y, method = "igd", ...)

Arguments

Value

a new statespace object with new values in kalman_params

Design a State-Space Model

Description

The function statespace builds a state-space model, with known or unknown variances. By default, this function builds a state-space model in the static setting, with a constant state (zero state noise covariance matrix) and constant observation noise variance.

Usage

statespace(X, y, kalman_params = NULL, viking_params = NULL, ...)

Arguments

Value

a statespace object.

Examples

set.seed(1)
### Simulate data
n <- 1000
d <- 5
Q <- diag(c(0,0,0.25,0.25,0.25))
sig <- 1

X <- cbind(matrix(rnorm((d-1)*n,sd=1),n,d-1),1)
theta <- matrix(rnorm(d), d, 1)
theta_arr <- matrix(0, n, d)
for (t in 1:n) {
  theta_arr[t,] <- theta
  theta <- theta + matrix(mvtnorm::rmvnorm(1,matrix(0,d,1),Q),d,1)
}
y <- rowSums(X * theta_arr) + rnorm(n, sd=sig)

####################
### Kalman Filter
# Default Static Setting
ssm <- viking::statespace(X, y)
plot(ssm)

# Known variances
ssm <- viking::statespace(X, y, kalman_params = list(Q=Q, sig=sig))
plot(ssm)

Viking: Variational bayesIan variance tracKING

Description

viking is the state-space estimation realized by Viking, generalizing the Kalman Filter to variance uncertainty.

Usage

viking(
  X,
  y,
  theta0,
  P0,
  hata0,
  s0,
  hatb0,
  Sigma0,
  n_iter = 2,
  mc = 10,
  rho_a = 0,
  rho_b = 0,
  learn_sigma = TRUE,
  learn_Q = TRUE,
  K = NULL,
  mode = "diagonal",
  thresh = TRUE,
  phi = logt,
  phi1 = logt1,
  phi2 = logt2
)

Arguments

Value

a list composed of the evolving value of all the parameters: theta_arr, P_arr, q_arr, hata_arr, s_arr, hatb_arr, Sigma_arr

References

J. de Vilmarest, O. Wintenberger (2021), Viking: Variational Bayesian Variance Tracking. <arXiv:2104.10777>