# Problem with definition variables with state space expectation?

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)

## Weird

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

In reply to Weird by mhunter

## Fixed

Log in or register to post comments

In reply to Fixed by mhunter

## Ok great :) But since 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

In reply to Ok great :) But since this by CharlesD

## Documentation fix

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

In reply to Documentation fix by mhunter

## Ok, good to know then, all

Log in or register to post comments