Last chance! 50% off unlimited learning
Sale ends in
mxExpectationStateSpaceContinuousTime(A, B, C, D, Q, R, x0, P0, u, t = NA,
dimnames = NA, thresholds = NA, threshnames = dimnames,
..., scores=FALSE)
mxExpectationSSCT(A, B, C, D, Q, R, x0, P0, u, t = NA,
dimnames = NA, thresholds = NA, threshnames = dimnames,
..., scores=FALSE)
mxExpectationStateSpaceContinuousTime
and mxExpectationSSCT
functions are identical. The latter is simply an abbreviated name. When using the former, tab completion is strongly encouraged to save tedious typing. Both of these functions are wrappers for the mxExpectationStateSpace function, which could be used for both discrete and continuous time modeling. However, there is a strong possibility of misunderstanding the model parameters when switching between discrete time and continuous time. The expectation matrices have the same names, but mean importantly different things so caution is warranted. The best practice is to use mxExpectationStateSpace
for discrete time models, and mxExpectationStateSpaceContinuousTime
for continuous time models. Expectation functions define the way that model expectations are calculated. That is to say, expectation functions define how a set of model matrices get turned into expectations for the data. When used in conjunction with the mxFitFunctionML, the mxExpectationStateSpace uses maximum likelihood prediction error decomposition (PED) to obtain estimates of free parameters in a model of the raw MxData object. Continuous time state space expectations treat the raw data as a multivariate time series of possibly unevenly spaced times with each row corresponding to a single occasion. Continuous time state space expectations implement a hybrid Kalman filter to produce expectations. The hybrid Kalman filter uses a Kalman-Bucy filter for the prediction step and the classical Kalman filter for the update step. It is a hybrid between the classical Kalman filter used for the discrete (but possibly unequally spaced) measurement occastions and the continous time Kalman-Bucy filter for latent variable predictions. Missing data handling is implemented in the same fashion as full information maximum likelihood for partially missing rows of data. Additionally, completely missing rows of data are handled by only using the prediction step from the Kalman-Bucy filter and omitting the update step. This model uses notation for the model matrices commonly found in engineering and control theory. The 'A', 'B', 'C', 'D', 'Q', 'R', 'x0', and 'P0' arguments must be the names of MxMatrix or MxAlgebraobjects with the associated properties of the A, B, C, D, Q, R, x0, and P0 matrices in the state space modeling approach. The 't' matrix must be a 1x1 matrix using definition variables that gives the times at which measurements occurred. The state space expectation is defined by the following model equations. #------------------------------------------------------------------------------
# Example 1
# Undamped linear oscillator, i.e. a noisy sine wave.
# Measurement error, but no dynamic error, single indicator.
# This example works great.
#--------------------------------------
# Data Generation
require(OpenMx)
set.seed(405)
tlen <- 200
t <- seq(1.2, 50, length.out=tlen)
freqParam <- .5
initialCond <- matrix(c(2.5, 0))
x <- initialCond[1,1]*cos(freqParam*t)
plot(t, x, type='l')
measVar <- 1.5
y <- cbind(obs=x+rnorm(tlen, sd=sqrt(measVar)), tim=t)
plot(t, y[,1], type='l')
#--------------------------------------
# Model Specification
#Note: the bounds are here only to keep SLSQP from
# stepping too far off a cliff. With the bounds in
# place, SLSQP finds the right solution. Without
# the bounds, SLSQP goes crazy.
cdim <- list('obs', c('ksi', 'ksiDot'))
amat <- mxMatrix('Full', 2, 2, c(FALSE, TRUE, FALSE, TRUE), c(0, -.1, 1, -.2),
name='A', lbound=-10)
bmat <- mxMatrix('Zero', 2, 1, name='B')
cmat <- mxMatrix('Full', 1, 2, FALSE, c(1, 0), name='C', dimnames=cdim)
dmat <- mxMatrix('Zero', 1, 1, name='D')
qmat <- mxMatrix('Zero', 2, 2, name='Q')
rmat <- mxMatrix('Diag', 1, 1, TRUE, .4, name='R', lbound=1e-6)
xmat <- mxMatrix('Full', 2, 1, TRUE, c(0, 0), name='x0', lbound=-10, ubound=10)
pmat <- mxMatrix('Diag', 2, 2, FALSE, 1, name='P0')
umat <- mxMatrix('Zero', 1, 1, name='u')
tmat <- mxMatrix('Full', 1, 1, name='time', labels='data.tim')
osc <- mxModel("LinearOscillator",
amat, bmat, cmat, dmat, qmat, rmat, xmat, pmat, umat, tmat,
mxExpectationSSCT('A', 'B', 'C', 'D', 'Q', 'R', 'x0', 'P0', 'u', 'time'),
mxFitFunctionML(),
mxData(y, 'raw'))
oscr <- mxRun(osc)
#--------------------------------------
# Results Examination
summary(oscr)
(ssFreqParam <- mxEval(sqrt(-A[2,1]), oscr))
freqParam
(ssMeasVar <- mxEval(R, oscr))
measVar
dampingParam <- 0
(ssDampingParam <- mxEval(-A[2,2], oscr))
dampingParam
Run the code above in your browser using DataLab