Chains don't move at all for AR(1) model -- please help!!

[CFQ]

Hi all,

I'm trying to fit a simple AR(1) model as below. But the chains don't move at all. According to the warning message, it suggests to increase adapt_delta. I've made it to 0.999 already but didn't solve the problem. Attached is the trace plot I got...

Please kindly help and feel free to let me know any suggestions that could make the chains move around. Thank you!!

Side note: I'm using

rstan (Version 2.8.0, packaged: 2015-09-19 14:48:38 UTC, GitRev: 05c3d0058b6a)

R version 3.1.2 (2014-10-31) -- "Pumpkin Helmet"

Platform: x86_64-apple-darwin13.4.0 (64-bit)

############ R code starts ############

library(rstan)

A <- 4

B <- 5

N <- A * B

a_n <- rep(1:A, each = B)

b_n <- rep(1:B, times = A)

nu <- 3

rho <- 0.9

sigma_eps <- 0.03

set.seed(1234)

logy <- rnorm(N, mean = 0, sd = 0.05)

y <- exp(logy)

logy_a <- logy[1:A]

X_lower <- log(0.5)

X_upper <- log(2)

test_data <- list(A = A, B = B, N = N,

                  y = y, logy = logy, logy_a = logy_a,

                  nu = nu, rho = rho, sigma_eps = sigma_eps,

                  a_n = a_n, b_n = b_n,

                  X_lower = X_lower, X_upper = X_upper)

test_model3 <- "

data {

int<lower=1> A;

real<lower=0, upper=1> rho;

real<lower=0> sigma_eps;

real logy_a[A];

real X_lower;

real X_upper;

}//end of data

parameters {

real<lower=X_lower, upper=X_upper> X_a[A];

}//end of parameters

model {

real sigmaSQ_eps_stat;

real sigma_eps_stat;

sigmaSQ_eps_stat <- sigma_eps^2 / (1 - rho^2);

sigma_eps_stat <- sqrt(sigmaSQ_eps_stat);

X_a[1] ~ normal(0, sigma_eps_stat) T[X_lower, X_upper];

for (a in 2:A) {

X_a[a] ~ normal(rho * X_a[a-1], sqrt(sigmaSQ_eps_stat * (1 - rho^2))) T[X_lower, X_upper];

}//end of a loop

// data model //

for (a in 1:A) {

logy_a[a] ~ normal(X_a[a], 0.5);

}//end of n loop

}//end of model

"

fit_test3 <- sampling(stan_model(model_code = test_model3),

                      data = test_data,

                      init = 0,

                      control = list(adapt_delta = 0.999,#default=0.8

                                     max_treedepth = 12),

                      iter = 2000, chains = 2, warmup = 2)

traceplot(fit_test3, par = dimnames(fit_test3)[[3]][sample(1:dim(fit_test3)[3], size = 5)])

############ R code ends ############

[Ben Goodrich]

On Tuesday, November 3, 2015 at 9:10:36 AM UTC-5, CFQ wrote:

Please kindly help and feel free to let me know any suggestions that could make the chains move around. Thank you!!

You basically have to Matt trick AR(1) processes for unknown parameters in Stan to get reasonable sampling. Like this:

parameters{
  real Z_a[A];
}
transformed parameters {
  real X_a[A];
  X_a[1] <- X_lower + inv_logit(sigma_eps * Z_a[1]) * (X_upper - X_lower);
  for (i in 2:A) {
    X_a[i] <- X_lower + inv_logit(rho * X_a[i-1] * sigma_eps * Z_a[i]) * (X_upper - X_lower);
  }
}
parameters {
  Z_a ~ normal(0,1);
  logy_a ~ normal(X_a, 0.5);
}

Ben

[Fengqing Chao]

Hi Ben,

Thank you so much for your suggestion! I just tried it using the updated model (attached). But I still need some help.

1. Z_a mix well, but X_a still don't mix. Please see attached trace plot.

2. In order to get the chain mix, I have to increase adapt_alpha instead of using the default. This is a very simple model with few parameters just to illustrate my problem. I was wandering if the number of parameter increases significantly, the chains won't mix well even if I set adapt_alpha very close to 1.

Thanks.

Regards,

CFQ

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/w6Tcdbar9zk/unsubscribe.
To unsubscribe from this group and all its topics, 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.

[Ben Goodrich]

On Tuesday, November 3, 2015 at 11:23:33 AM UTC-5, CFQ wrote:

Thank you so much for your suggestion! I just tried it using the updated model (attached). But I still need some help.

1. Z_a mix well, but X_a still don't mix. Please see attached trace plot.

2. In order to get the chain mix, I have to increase adapt_alpha instead of using the default. This is a very simple model with few parameters just to illustrate my problem. I was wandering if the number of parameter increases significantly, the chains won't mix well even if I set adapt_alpha very close to 1.

If you have to set adapt_delta very close to 1 to prevent divergent transitions, then so be it no matter what the mixing rate is. But in your R code, logy_a is not generated as a function of X_a, so there is no information in the data with which to estimate it.

Ben

[Andrew Gelman]

Worth putting this in the manual, I think.

A

--
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.

[andre.p...@googlemail.com]

This model confuses me. 

In R you write:

logy_a <- rnorm(A, mean = 0, sd = 0.05)

In STAN you write:

logy_a ~ normal(X_a, 0.5);

Why there are different sds?

Further I don't see an AR(1) process in R code.

Andre

[Michael Betancourt]

Definitely — Ben have you seen a trade off in performance between

this and the standard parameterization when the data becomes more

dense?  I can fit simple models using the standard parameterization

no problem but I have definitely run into problems with larger models.

[Fengqing Chao]

Hi Andre,

logy_a <- rnorm(A, mean = 0, sd = 0.05)

The line above is to generate the input data for STAN model.

logy_a ~ normal(X_a, 0.5);

The line above in STAN is the data model, where I assume that logy_a follows an AR(1). X_a is the parameter vector that I want to estimate. As you mentioned, yes, I didn't generate logy_a by AR(1). The AR(1) process is coded in STAN model code. What I want to do is to make logy_a resemble some real-life data on log scale. I don't see the necessity of generating input data logy_a in R using AR(1) process. Whatever the underlying distribution it has, the chains should at least move around instead of staying at the exact same value.

Regards,

CFQ

[andre.p...@googlemail.com]

Hi CFQ,

 I think your parameter array (vector) is Z_a, since its defined as parameter in the parameter section and

X_a is an input. Sorry, but an AR(1) model normal has constant, and a rho, phi, or whatever you call that correlates to the past X_t-1

and if you want a sigma. As I see, you defined A parameters, since A=400, its 400+ parameters. Did you have had a look in the manual. 

There is an AR(1) model example, that is working well. Why should the chain move? Isn't it you defined something like a table lookup? 

Sorry, I really don't want to sound harsh. I just get tons of questions in my brain, and maybe you trying something genius, but 

I cannot follow. 

Best regards,

Andre

[CFQ]

Hi Ben,

I re-write the STAN code as below (I edit a little bit about the inv_logit transformation you suggested earlier), and now both the chains for X_a and Z_a move around!

############ R code starts ############

library(rstan)

A <- 400

rho <- 0.9
sigma_eps <- 0.03
set.seed(1234)

logy_a <- rnorm(A, mean = 0, sd = 0.05)
y_a <- exp(logy_a)

X_lower <- log(0.5)
X_upper <- log(2)

test_data <- list(A = A, logy_a = logy_a,

                  rho = rho, sigma_eps = sigma_eps,

                  X_lower = X_lower, X_upper = X_upper)
test_model3 <-

"

data {
int<lower=1> A;
real<lower=0, upper=1> rho;
real<lower=0> sigma_eps;
real logy_a[A];
real X_lower;
real X_upper;
}//end of data

parameters{
  real Z_a[A];
}
transformed parameters {

  real<lower=X_lower, upper=X_upper> X_a[A];
real<lower=0> sigma_eps_stat;


sigma_eps_stat <- sqrt(sigma_eps^2 / (1 - rho^2));


  X_a[1] <- X_lower + inv_logit(sigma_eps_stat * Z_a[1]) * (X_upper - X_lower);

  for (i in 2:A) {
    X_a[i] <- X_lower +

              inv_logit(rho * X_a[i-1] + sigma_eps * Z_a[i]) * (X_upper - X_lower);
  }
}
model {

  Z_a ~ normal(0,1);
  logy_a ~ normal(X_a, 0.5);
}
"

fit_test3 <- sampling(stan_model(model_code = test_model3),
                      data = test_data,
                      init = 0,
                      control = list(adapt_delta = 0.999,#default=0.8
                                     max_treedepth = 12),

                      iter = 1000, chains = 2, warmup = 2)
traceplot(fit_test3, par = c(paste0("X_a[", c(10, 211, 352), "]"),
                             paste0("Z_a[", c(10, 211, 352), "]")))

############ R code ends ############

[CFQ]

Hi Andre,

Thanks for the discussion! Here're my thoughts:

I think your parameter array (vector) is Z_a, since its defined as parameter in the parameter section and
X_a is an input.

##cfq: both Z_a and X_a are parameters. X_a is in the "transformed parameter" block. And X_a is my parameter of interest. Z_a is the inverse logit transformation of X_a as Ben suggested in order to get reasonable sampling.

Sorry, but an AR(1) model normal has constant, and a rho, phi, or whatever you call that correlates to the past X_t-1
and if you want a sigma.

##cfq: I agree and see nothing wrong in my model regarding this point.

you defined A parameters, since A=400, its 400+ parameters.

##cfq: compare to my actual model, the number of parameter in this example is quite small.

 Did you have had a look in the manual. 
There is an AR(1) model example, that is working well.

##cfq: yeah thanks. But what I need are not only the rho, sigma, but also the time-series vector itself, X_a. That's not necessary in my example model, but I do need to save this X_a parameter in my actual model which is a time-series plus something else.

Why should the chain move? Isn't it you defined something like a table lookup? 

 ##cfq: The chain moves so that it can explore the whole range of prior and get reasonable posterior distribution. unless I purposely model the parameter in a way that it's like a constant but I didn't. So I'm not sure what you mean by "a table lookup".

CFQ

[Ben Goodrich]

On Wednesday, November 4, 2015 at 8:01:10 AM UTC-5, Michael Betancourt wrote:

Definitely — Ben have you seen a trade off in performance between

this and the standard parameterization when the data becomes more

dense?  I can fit simple models using the standard parameterization

no problem but I have definitely run into problems with larger models.

My memory is consistent with that but I haven't done a AR(1) model in a while. I think I was doing autoregressive coefficients with the General Social Survey data.

Ben
 

On Nov 3, 2015, at 10:52 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

Worth putting this in the manual, I think.

A

On Nov 3, 2015, at 10:00 AM, Ben Goodrich <goodri...@gmail.com> wrote:

On Tuesday, November 3, 2015 at 9:10:36 AM UTC-5, CFQ wrote:

Please kindly help and feel free to let me know any suggestions that could make the chains move around. Thank you!!

You basically have to Matt trick AR(1) processes for unknown parameters in Stan to get reasonable sampling. Like this:

parameters{
  real Z_a[A];
}
transformed parameters {
  real X_a[A];
  X_a[1] <- X_lower + inv_logit(sigma_eps * Z_a[1]) * (X_upper - X_lower);
  for (i in 2:A) {
    X_a[i] <- X_lower + inv_logit(rho * X_a[i-1] * sigma_eps * Z_a[i]) * (X_upper - X_lower);
  }
}
parameters {
  Z_a ~ normal(0,1);
  logy_a ~ normal(X_a, 0.5);
}


Ben

--
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+unsubscribe@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+unsubscribe@googlegroups.com.