Login

News

Problem with definition variables with state space expectation?

Posted on
No user picture. CharlesD Joined: 04/30/2013
Forums
OpenMx General Help
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)

Replied on Mon, 02/02/2015 - 11:53
Picture of user. mhunter Joined: 07/31/2009

What you're writing should work, but I'm replicating the same behavior. When you use the pre-run model, these evaluations work.


mxEval(data.V1, ssModel, compute=T, defvar.row=1)
mxEval(u, ssModel, compute=T, defvar.row=1)

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!

Log in or register to post comments
Replied on Mon, 02/02/2015 - 19:20
Picture of user. mhunter Joined: 07/31/2009

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!
Log in or register to post comments
Replied on Tue, 02/03/2015 - 05:31
No user picture. CharlesD Joined: 04/30/2013

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]<-A*x[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)

Log in or register to post comments
Replied on Tue, 02/03/2015 - 10:51
Picture of user. mhunter Joined: 07/31/2009

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.

Log in or register to post comments
  Login

News