Help with rstimation of dlm with unknown parameters in transition equation

[Peyman Noursalehi]

Hi,

I'm reading the Dynamic Linear Models with R book, where most of chapter 4 is devoted to Bayesian estimation of parameters. They code most of it manually though, and it can get quite tricky for complicated models. I'm trying to reproduce the example on p.186 in Stan, where the model to be estimated consists of local linear trend and AR(2) components

Priors for (ϕ1ϕ2) are N(0,(2/3)^2) and N(0,(1/3)^2) and the inverses of variances are assumed to be independent with Gamma priors g(a2/b,a/b) where a = 1 and b =1000.

I found some discussions here before on constraining parameters to be stationary, so that part has been partially resolved. 

However, the estimation fails because "error occurred during calling the sampler; sampling not done" after a ton of messages that complain "Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: validate transformed params: W is not positive definite". 

I'm a beginner in Stan and I'm certain that I've messed up somewhere in the code. Any help is appreciated!

Thanks,

Peyman

[Jeffrey Arnold]

The reason this isn't working is that you set the type of W to be a covariance matrix, and W is not positive definite (it has a 0 on its diagonal). It should work if you replace "cov_matrix [4] W"; with "matrix [4, 4] W;"

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Michael Betancourt]

That will still fail when you try to evaluate the multivariate normal density for Omega_{t}.

Firstly you don’t want to treat W is a full covariance matrix given that it’s just a diagonal,

instead decompose it into a bunch of independent one-dimensional normal distributions.

The 0 variance for the last element of Omega_{t} implies that it is a constant, so you want

to pull that out of the parameters entirely, leaving just a 3-dimensional vector of Omega_{t}

that do have a well-posed distribution.

[Peyman]

Thanks for your suggestion. I first changed the type to matrix. It just never converges "Warning messages:1: There were 3970 divergent transitions after warmup." I've tried changing priors, but the probelms "Exception thrown at line 90: gaussian_dlm_obs_log: G[11] is -1.#INF, but must be finite! " and "Exception thrown at line 90: gaussian_dlm_obs_log: V[1] is 1.#INF, but must be finite! " happen show up a lot in the log. 

data {

  int <lower = 0> N;

  matrix [1, N] y;

}

transformed data {

  matrix [4, 1] F;

  vector [4] m0;

  cov_matrix [4] C0;

  

  F [1, 1] <- 1;

  F [2, 1] <- 0;

  F [3, 1] <- 1;

  F [4, 1] <- 0;

  m0 [1] <- 0;

  m0 [2] <- 0;

  m0 [3] <- 0;

  m0 [4] <- 0;

  

  C0 <- diag_matrix(rep_vector(1.0e+7, 4));

}

parameters {

  // need to impose stationarity constraints

 real <lower=-1,upper=1> phi2; 

 real <upper=(1 - fabs(phi2))> phi1; 

 

 

  real <lower = 0> sigma ; // for V

  // real<lower = 0> W1;

  // real<lower = 0> W2;

  // real<lower = 0> W3;

  

  vector<lower = 0>[4] W_diag;

  

}

transformed parameters {

  vector [1] V;

  matrix [4, 4] W;

  matrix [4, 4] G;

  G [1, 1] <- 1;

  G [1, 2] <- 1;

  G [1, 3] <- 0;

  G [1, 4] <- 0;

  G [2, 1] <- 0;

  G [2, 2] <- 1;

  G [2, 3] <- 0;

  G [2, 4] <- 0;  

  G [3, 1] <- 0;

  G [3, 2] <- 0;

  G [3, 3] <- phi1;

  G [3, 4] <- 1;

  G [4, 1] <- 0;

  G [4, 2] <- 0;  

  G [4, 3] <- phi2;

  G [4, 4] <- 0;

  

  

  V [1] <- sigma  * sigma ;

  W <- diag_matrix(W_diag);

  // W [1, 1] <- W1;

  // W [1, 2] <- 0;

  // W [1, 3] <- 0;

  // W [1, 4] <- 0;

  // W [2, 1] <- 0;

  // W [2, 2] <- W2;

  // W [2, 3] <- 0;

  // W [2, 4] <- 0;  

  // W [3, 1] <- 0;

  // W [3, 2] <- 0;

  // W [3, 3] <- W3;

  // W [3, 4] <- 1;

  // W [4, 1] <- 0;

  // W [4, 2] <- 0;  

  // W [4, 3] <- 0;

  W [4, 4] <- 0;

}

model {

  sigma ~ uniform (0, 5);

  phi2 ~ normal(0,2.0/3);

  phi1 ~ normal(0,1.0/3);

  W[1,1] ~ inv_gamma(0.001, 0.001);

  W[2,2] ~ inv_gamma(0.001, 0.001);

  W[3,3] ~ inv_gamma(0.001, 0.001);

  // W1 ~ inv_gamma(1, 10);

  // W2 ~ inv_gamma(1, 10);

  // W3 ~ inv_gamma(1, 10);

  y ~ gaussian_dlm_obs (F, G, V, W, m0, C0);

  

}

[Peyman]

I think I have changed my code to do it as you suggested:
But it quickly quits saying "Stan model 'linear_trend_ar2_model' does not contain samples." which I believe means I've misspecified my model

data {

  int <lower = 0> N;

  matrix [1, N] y;

}

transformed data {

  matrix [4, 1] F;

  vector [4] m0;

  cov_matrix [4] C0;

  

  F [1, 1] <- 1;

  F [2, 1] <- 0;

  F [3, 1] <- 1;

  F [4, 1] <- 0;

  m0 [1] <- 0;

  m0 [2] <- 0;

  m0 [3] <- 0;

  m0 [4] <- 0;

  

  C0 <- diag_matrix(rep_vector(1.0e+7, 4));

}

parameters {

  // need to impose stationarity constraints

 real <lower=-1,upper=1> phi2; 

 real <upper=(1 - fabs(phi2))> phi1; 

 

 

  real <lower = 0> sigma ; // for V

  real<lower = 0> W1;

  real<lower = 0> W2;

  real<lower = 0> W3;

  

  // vector<lower = 0>[4] W_diag;

  

}

transformed parameters {

  vector [1] V;

  matrix [4, 4] W;

  matrix [4, 4] G;

  G [1, 1] <- 1;

  G [1, 2] <- 1;

  G [1, 3] <- 0;

  G [1, 4] <- 0;

  G [2, 1] <- 0;

  G [2, 2] <- 1;

  G [2, 3] <- 0;

  G [2, 4] <- 0;  

  G [3, 1] <- 0;

  G [3, 2] <- 0;

  G [3, 3] <- phi1;

  G [3, 4] <- 1;

  G [4, 1] <- 0;

  G [4, 2] <- 0;  

  G [4, 3] <- phi2;

  G [4, 4] <- 0;

  

  

  V [1] <- sigma  * sigma ;

  // W <- diag_matrix(W_diag);

  W [1, 1] <- W1;

  W [1, 2] <- 0;

  W [1, 3] <- 0;

  W [1, 4] <- 0;

  W [2, 1] <- 0;

  W [2, 2] <- W2;

  W [2, 3] <- 0;

  W [2, 4] <- 0;

  W [3, 1] <- 0;

  W [3, 2] <- 0;

  W [3, 3] <- W3;

  W [3, 4] <- 1;

  W [4, 1] <- 0;

  W [4, 2] <- 0;

  W [4, 3] <- 0;

  W [4, 4] <- 0;

}

model {

  sigma ~ uniform (0, 5);

  phi2 ~ normal(0,2.0/3);

  phi1 ~ normal(0,1.0/3);

  // W[1,1] ~ inv_gamma(0.001, 0.001);

  // W[2,2] ~ inv_gamma(0.001, 0.001);

  // W[3,3] ~ inv_gamma(0.001, 0.001);

  W1 ~ inv_gamma(1, 0.001);

  W2 ~ inv_gamma(1, 0.001);

  W3 ~ inv_gamma(1, 0.001);

  y ~ gaussian_dlm_obs (F, G, V, W, m0, C0);

  

}

[Bob Carpenter]

The problem may be that you're trying to give W an
all-zero row.

You want an upper bound on sigma to match your uniform
distribution (another one our more verbose warnings would
have caught). But that's probably not your issue.

Did you want upper and lower bounds on the phi?

You almost never want a 1 x N matrix in Stan --- y should
probably be a row vector.

Starting a covariance matrix here:

> C0 <- diag_matrix(rep_vector(1.0e+7, 4));

is not a good idea unless you're going to be
drawing values in the tens of millions at which point,
you probably want to think about rescaling somehow.

Also, it's almost never a good idea to use unit covariance
matrices in Stan because it winds up doing a lot of expensive
and unnecessary matrix algebra. But then we don't have
a specialized unit-diagonal form of that density.

- Bob

>> Priors for (ϕ1ϕ2) are N(0,(2/3)^2) and N(0,(1/3)^2) and the inverses of variances are assumed to be independent with Gamma priors g(a2/b,a/b) where a = 1 and b =1000.
>> I found some discussions here before on constraining parameters to be stationary, so that part has been partially resolved.
>> However, the estimation fails because "error occurred during calling the sampler; sampling not done" after a ton of messages that complain "Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: validate transformed params: W is not positive definite".
>>
>> I'm a beginner in Stan and I'm certain that I've messed up somewhere in the code. Any help is appreciated!
>>
>> Thanks,
>> Peyman
>>
>> --
>> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
>> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
>> To post to this group, send email to stan-...@googlegroups.com.
>> For more options, visit https://groups.google.com/d/optout.
>>
>>
>> --
>> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
>> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
>> To post to this group, send email to stan-...@googlegroups.com.
>> For more options, visit https://groups.google.com/d/optout.
>
>
> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

> <gdp5004.dat><linear_trend_ar2_model.stan><linear_trend_ar2_model_for_stan.r>