Estep: expectation step in sparse expectation-maximization algorithm

In hdiVAR: Statistical Inference for Noisy Vector Autoregression

expectation step in sparse expectation-maximization algorithm

Description

Compute conditional expectation and covariance of x_t given y_1,\ldots,y_T and current parameter estimates of A, \sigma_\eta,\sigma_\epsilon via kalman filter and smoothing.

Usage

Estep(Y,A_init,sig_eta_init,sig_epsilon_init,X_init,P_init)

Arguments

Y	observations of time series, a p by T matrix.
A_init	current estimate of transition matrix A.
sig_eta_init	current estiamte of \sigma_\eta.
sig_epsilon_init	current estiamte \sigma_\epsilon.
X_init	current estimate of latent x_1 at the first iteration, a p-dimensional vector.
P_init	current covariance estimate of latent x_1 at the first iteration, a p by p matrix.

Value

a list of conditional expectations and covariances for the sequential Maximization step.

EXtT	a p by T matrix of column E[x_t | y_1,\ldots,y_T, \hat{A},\hat{\sigma}_\eta,\hat{\sigma}_\epsilon].
EXtt	a p by p by T tensor of first-two-mode slice E[x_t x_t^\top | y_1,\ldots,y_T, \hat{A},\hat{\sigma}_\eta,\hat{\sigma}_\epsilon].
EXtt1	a p by p by T-1 matrix of first-two-mode slice E[x_t x_{t+1}^\top | y_1,\ldots,y_T, \hat{A},\hat{\sigma}_\eta,\hat{\sigma}_\epsilon].

Author(s)

Xiang Lyu, Jian Kang, Lexin Li

Examples

p= 2; Ti=10  # dimension and time
A=diag(1,p) # transition matrix
sig_eta=sig_epsilon=0.2 # error std
Y=array(0,dim=c(p,Ti)) #observation t=1, ...., Ti
X=array(0,dim=c(p,Ti)) #latent t=1, ...., T
Ti_burnin=100 # time for burn-in to stationarity
for (t in 1:(Ti+Ti_burnin)) {
  if (t==1){
    x1=rnorm(p)
  } else if (t<=Ti_burnin) { # burn in
    x1=A%*%x1+rnorm(p,mean=0,sd=sig_eta)
  } else if (t==(Ti_burnin+1)){ # time series used for learning
    X[,t-Ti_burnin]=x1
    Y[,t-Ti_burnin]=X[,t-Ti_burnin]+rnorm(p,mean=0,sd=sig_epsilon)
  } else {
    X[,t- Ti_burnin]=A%*%X[,t-1- Ti_burnin]+rnorm(p,mean=0,sd=sig_eta)
    Y[,t- Ti_burnin]=X[,t- Ti_burnin]+rnorm(p,mean=0,sd=sig_epsilon)
  }
}

Efit=Estep(Y,A,sig_eta,sig_epsilon,x1,diag(1,p))

hdiVAR documentation built on May 31, 2023, 7:27 p.m.