Lets say I have the following state space model:

$y_t = \beta_t x_t + \epsilon_t$

$\beta_{t+1} = \mu_t + \beta_t \eta_t$

$\mu_{t+1} = \mu_t + \omega_t$

All my true values for $y$ are known, but I want to obtain fitted values. In Jags I run the following.

# The data (an any other input) we pass to jags
jags.data = list("Y"=y, "N"=N, Y1=y[1], C=x1)
# The parameters that we are monitoring (must monitor at least 1)
jags.params=c("sd.q","sd.r","Y", "X1", "X2")
mod_ss = jags(jags.data, parameters.to.save=jags.params, model.file=model.loc, n.chains = 3, 
              n.burnin=5000, n.thin=1, n.iter=10000, DIC=TRUE)  

attach.jags(mod_ss, overwrite=TRUE);

But then the values that I get back for $y$ are the same as the true values. I sort of get this, because Jags can of course not find better values than the true values and its Bayesian, but is there not a way to let it fit new values? Such that my true $y$ values are used for estimation but won't be the same as the estimated $y$?

Edit: Jags code

# jags model specification
jagsscript = cat("
  model {  
    # priors on parameters
    # Make sure mu prior is scaled to the data
    mu ~ dnorm(Y1, 1/(Y1*100)); 
  tau.q ~ dgamma(0.001,0.001); # This is inverse gamma
    sd.q <- 1/sqrt(tau.q); # sd is treated as derived parameter
  tau.r ~ dgamma(0.001,0.001); # This is inverse gamma
    sd.r <- 1/sqrt(tau.r); # sd is treated as derived parameter
    u ~ dnorm(0, 0.01);

    # If X[0] = mu instead of X[1]
    # X1[1] ~ dnorm(mu+u, tau.q)
    X1[1] <- mu;
  X2[1] <- mu;
    Y[1] ~ dnorm(X1[1]*C[1], tau.r);
    # Jags is not vectorized, so we have to loop over observations
    for(i in 2:N) {
    predX1[i] <- X1[i-1]+u; 
  predX2[i] <- X2[i-1]+u; 
    X1[i] ~ dnorm(predX1[i] + predX2[i], tau.q); # Process variation
  X2[i] ~ dnorm(predX2[i], tau.q); # Process variation
    Y[i]  ~ dnorm(X1[i]*C[i], tau.r); # Observation variation

  • $\begingroup$ What is your JAGS code? $\endgroup$
    – Tim
    Apr 12 '15 at 10:16
  • $\begingroup$ @Tim Please see the edit that I just made. Thanks. $\endgroup$
    – student07
    Apr 12 '15 at 10:50

Add a lines along the lines of:

Y_pred[1] ~ dnorm(X1[1]*C[1], tau.r);

and further down:

Y_pred[i]  ~ dnorm(X1[i]*C[i], tau.r);

But don't "define" Y_pred (by adding it as data). Then Y_pred will contain draws from the predictive distribution of Y.

  • $\begingroup$ Ah, this is clever! I will try it out right away and see if it gives the desired output! Thank you very much! $\endgroup$
    – student07
    Apr 12 '15 at 11:47
  • $\begingroup$ It works now in the sense that I get values that are different from the true values. However, I mainly have two questions: 1. Y_pred now gives a matrix which is (13 x #number of obs), why is this so? I took the first row now but I am not sure which row I should take. 2. If I take the first row the values are nowhere near the actual values, i.e. the scale of the fitted values is nowhere far from correct, how could this be so? Thanks again! $\endgroup$
    – student07
    Apr 12 '15 at 12:13
  • $\begingroup$ Since Y depends on C, you get one Y[.] for every C[.]. Then Y_pred[1] is the prediction given C[1]. Shat do you mean with that the scale is different? If Y_pred is very spread out it just meand that there is very little information in the model+data and that Y_pred is uncertain. $\endgroup$ Apr 12 '15 at 14:38
  • $\begingroup$ Yes, but why do I get 13 rows? $\endgroup$
    – student07
    Apr 12 '15 at 14:56
  • $\begingroup$ You have 13 C, right? $\endgroup$ Apr 12 '15 at 15:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.