expectation_maximization {viking}

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

X

explanatory variables

y

time series

n_iter

number of iterations of the EM algorithm

Q_init

initial covariance matrix on the state noise

sig_init

(optional, default 1) initial value of the standard deviation of the observation noise

verbose

(optional, default 1000) frequency for prints

lambda

(optional, default 10^-9) regularization parameter to avoid singularity

mode_diag

(optional, default FALSE) if TRUE then we restrict the search to diagonal matrices for Q

p1

(optional, default 0) deterministic value of P1 = p1 I

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)