Am I doing something wrong? I was surprised when my model seemed to be ignoring my definition variables, but then I also notice that mxEval interprets them all as NA... this is using a build from yesterday, 1st Feb 2015, on 64 bit windows.
data(demoOneFactor)
nvar <- ncol(demoOneFactor)
varnames <- colnames(demoOneFactor)
demoOneFactorInputs <- cbind(demoOneFactor, V1=rnorm(nrow(demoOneFactor)))
ssModel <- mxModel(model="State Space Inputs Manual Example",
mxMatrix("Full", 1, 1, TRUE, .3, name="A"),
mxMatrix("Full", 1, 1, TRUE, values=1, name="B"),
mxMatrix("Full", nvar, 1, TRUE, .6, name="C", dimnames=list(varnames, "F1")),
mxMatrix("Zero", nvar, 1, name="D"),
mxMatrix("Diag", 1, 1, FALSE, 1, name="Q"),
mxMatrix("Diag", nvar, nvar, TRUE, .2, name="R"),
mxMatrix("Zero", 1, 1, name="x0"),
mxMatrix("Diag", 1, 1, FALSE, 1, name="P0"),
mxMatrix("Full", 1, 1, FALSE, labels="data.V1", name="u"),
mxData(observed=demoOneFactorInputs, type="raw"),
mxExpectationStateSpace("A", "B", "C", "D", "Q", "R", "x0", "P0", u="u"),
mxFitFunctionML()
)
ssRun <- mxRun(ssModel)
mxEval(data.V1, ssRun, compute=T, defvar.row=1)
Printer-friendly version
What you're writing should work, but I'm replicating the same behavior. When you use the pre-run model, these evaluations work.
And when you generate data based on a model with non-zero weights in the "B" and/or "D" matrices, the expectation correctly estimates the parameters. A model like this is in our test suite. http://openmx.psyc.virginia.edu/svn/trunk/models/passing/StateSpaceInputs.R
So, I'm very confident that internally the definition variables/inputs are being used correctly. There does seems to be a problem in evaluating them with mxEval after a model has been run.
I've created a Issue for this bug and I'll be fixing it soon. http://openmx.psyc.virginia.edu/issue/2015/02/mxeval-run-state-space-model
Thanks for finding it!
This ended up being a data sorting issue. It's now fixed in the latest svn (r4208) and will be part of the next binary release. Thanks again for finding this!
Ok great :) But since this was not the problem I was having, I worked out what was. According to my understanding, although you say that this:
x[t+1] = A x[t] + B u[t] + q[t]
is the state equation, it 'seems' to be functioning as this:
x[t+1] = A x[t] + B u[t+1] + q[t]
The following code hopefully shows either a) what I mean, or b) where my thinking is wrong!
A<-diag(.5,1)
B<-matrix(c(0,2.7),nrow=1)
Q<-diag(.1,1)
Tpoints<-200
u<-matrix(c(rep(0,Tpoints),sample(c(0,0,0,0,0,0,10),size=Tpoints,replace=T)),byrow=T,nrow=2)
q<-rnorm(200,0,sqrt(2))
x<-rep(NA,200)
x[1]<-1.5
for(i in 2:Tpoints){
x[i]<-Ax[i-1] + B%%u[,i-1,drop=F] + q[i-1]
}
plot(x,type='l')
data<-cbind(x,u[2,])
colnames(data)<-c('y','input')
model<-mxModel(
mxData(data,type='raw'),
mxMatrix(name='A',free=T,nrow=1,ncol=1),
mxMatrix(name='B',free=T,nrow=1,ncol=2),
mxMatrix(name='C',values=1,free=F,nrow=1,ncol=1,dimnames=list('y','x')),
mxMatrix(name='D',values=0,free=F,nrow=1,ncol=2),
mxMatrix(name='Q',free=T,nrow=1,ncol=1),
mxMatrix(name='x0',free=T,nrow=1,ncol=1),
mxMatrix(name='P0',values=1,free=F,nrow=1,ncol=1),
mxMatrix(name='Q',free=T,nrow=1,ncol=1),
mxMatrix(name='R',values=0,free=F,nrow=1,ncol=1),
mxMatrix(name='u',values=c(1,0),labels=c(NA,'data.input'),free=F,nrow=2,ncol=1),
mxExpectationStateSpace('A', 'B', 'C', 'D', 'Q', 'R', 'x0', 'P0', 'u', dimnames='y'),
mxFitFunctionML()
)
fit<-mxRun(model)
summary(fit)
now shift the u matrix in time
data[2:(Tpoints),'input']<-data[1:(Tpoints-1),'input']
model<-mxModel(
mxData(data,type='raw'),
mxMatrix(name='A',free=T,nrow=1,ncol=1),
mxMatrix(name='B',free=T,nrow=1,ncol=2),
mxMatrix(name='C',values=1,free=F,nrow=1,ncol=1,dimnames=list('y','x')),
mxMatrix(name='D',values=0,free=F,nrow=1,ncol=2),
mxMatrix(name='Q',free=T,nrow=1,ncol=1),
mxMatrix(name='x0',free=T,nrow=1,ncol=1),
mxMatrix(name='P0',values=1,free=F,nrow=1,ncol=1),
mxMatrix(name='Q',free=T,nrow=1,ncol=1),
mxMatrix(name='R',values=0,free=F,nrow=1,ncol=1),
mxMatrix(name='u',values=c(1,0),labels=c(NA,'data.input'),free=F,nrow=2,ncol=1),
mxExpectationStateSpace('A', 'B', 'C', 'D', 'Q', 'R', 'x0', 'P0', 'u', dimnames='y'),
mxFitFunctionML()
)
fit<-mxRun(model)
summary(fit)
You're right. The state equation that's used has the "instantaneous" input effects for the dynamics (and the measurement).
x[t] = A x[t-1] + B u[t] + q[t]
not the lagged input effects.
x[t+1] = A x[t] + B u[t] + q[t]
If you want the lagged input effects you'll have to create a lagged variable yourself. I'll update the documentation to reflect the true state of things.
Ok, good to know then, all working fine now.