Divergent transitions & BFMI low in a state-space model

ytan · August 13, 2017, 7:28am

I ran the following state space model, but got a warning message. Please find an attached file to reproduce my results.

stan code

data {
  int P;
  int T;
  int T_pred;
  real<lower=0, upper=10> H[P,T];
}

parameters {
  real mu_H[P,T];
  real<lower=0> s_mu[P];
  real<lower=0> s_H;
}

model {
  for (p in 1:P){
    s_mu[p] ~ student_t(4,0,2); 
  }
  for (t in 3:T)
    for (p in 1:P){
      mu_H[p,t] ~ normal(2 * mu_H[p,t - 1] - mu_H[p,t - 2], s_mu[p]);
      //assuming that s_mu differs across prefectures 
      H[p,t] ~ normal(mu_H[p,t], s_H);
    }
}

generated quantities {
  real mu_all[P,T + T_pred];
  real h_pred[P,T_pred];
  for (t in 1:T)
    for (p in 1:P)
      mu_all[p,t] = mu_H[p,t];
  for (t in 1:T_pred)
    for (p in 1:P){
      mu_all[p,T + t] = normal_rng(2 * mu_all[p,T + t - 1]
                      - mu_all[p,T + t - 2], s_mu[p]);
      h_pred[p,t] = normal_rng(mu_all[p,T + t], s_H);
    }
}

r code

setwd("~/docdis")
library(rstan)
library(dplyr)
library(tidyr)
#import data
d <- read.csv(file ='input/ch3/reginfotest.csv', header = TRUE,
                       stringsAsFactors = FALSE)

#long -> wide
d <- spread(d, key = year, value = mean)

data <- list(P = nrow(d),T = ncol(d) - 1,T_pred = 3,H = d[,-1])
stanmodel <- stan_model(file='model/ch3/model-treghapSNC_simpmatxp.stan')
fit <- sampling(
  stanmodel,
  data = data,
  seed =243,
  chains=4, iter=2000, warmup=1000, thin=1, control=list(adapt_delta=0.99, max_treedepth=15)
)

warning message

1: There were 68 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
Runtime warnings and convergence problems
2: There were 2 chains where the estimated Bayesian Fraction of Missing Information was low. See
Runtime warnings and convergence problems
3: Examine the pairs() plot to diagnose sampling problems

Then, I ran another stan code with reparameterization to solve them as follows, but got another warning message. I suppose that this message is related to mistakenly specified index, but I could not find any solutions for this. I would appreciate if someone would respond to me.

stan code

 data {
   int P;
   int T;
   int T_pred;
   real<lower=0, upper=10> H[P,T];
 }
 
 parameters {
   real<lower=0> s_mu;
   real mu_H_raw[P,T-2];
   real<lower=0> s_H;
 }
 
 transformed parameters {
   real mu_H[P,T];
   for (t in 3:T)
     for (p in 1:P)
       mu_H[p,t] = 2 * mu_H[p,t - 1] - mu_H[p,t - 2]
                 + s_mu * mu_H_raw[p,t]; //reparameterization
 }
 model {
   for (t in 3:T)
     for (p in 1:P){
       mu_H_raw[p,t] ~ normal(0, 1); //reparameterization
       H[p,t] ~ normal(mu_H[p,t], s_H);
     }
 }
 
 generated quantities {
   real mu_all[P,T + T_pred];
   real h_pred[P,T_pred];
   for (t in 1:T)
     for (p in 1:P)
       mu_all[p,t] = mu_H[p,t];
   for (t in 1:T_pred)
     for (p in 1:P){
       mu_all[p,T + t] = 2 * mu_all[p,t - 1] - mu_all[p,t - 2]
                         + s_mu * normal_rng(0, 1);
       h_pred[p,t] = normal_rng(mu_all[p,T + t], s_H);
     }
 }

r code

setwd("~/docdis")
library(rstan)
library(dplyr)
library(tidyr)
#import data
d <- read.csv(file ='input/ch3/reginfotest.csv', header = TRUE,
                       stringsAsFactors = FALSE)

#long -> wide
d <- spread(d, key = year, value = mean)

data <- list(P = nrow(d),T = ncol(d) - 1,T_pred = 3,H = d[,-1])
stanmodel <- stan_model(file='model/ch3/model-treghapSNC_simpmatxrep.stan')
fit <- sampling(
  stanmodel,
  data = data,
  seed = 243,
  chains=4, iter=2000, warmup=1000, thin=1
)

warning message

SAMPLING FOR MODEL ‘model-treghapSNC_simpmatxrep’ NOW (CHAIN 1).
Unrecoverable error evaluating the log probability at the initial value.
Exception: : accessing element out of range. index 8 out of range; expecting index to be between 1 and 7; index position = 2mu_H_raw (in ‘model250c34c853fe_model_treghapSNC_simpmatxrep’ at line 18)

[1] “Error in sampler$call_sampler(args_list[[i]]) : "
[2] " Exception: : accessing element out of range. index 8 out of range; expecting index to be between 1 and 7; index position = 2mu_H_raw (in ‘model250c34c853fe_model_treghapSNC_simpmatxrep’ at line 18)”
[1] “error occurred during calling the sampler; sampling not done”

session info

R version 3.4.1 (2017-06-30)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

locale:
[1] LC_COLLATE=Japanese_Japan.932 LC_CTYPE=Japanese_Japan.932
[3] LC_MONETARY=Japanese_Japan.932 LC_NUMERIC=C
[5] LC_TIME=Japanese_Japan.932

attached base packages:
[1] stats graphics grDevices utils datasets methods base

other attached packages:
[1] ggmcmc_1.1 tidyr_0.6.3 dplyr_0.7.1
[4] rstan_2.16.2 StanHeaders_2.16.0-1 ggplot2_2.2.1

loaded via a namespace (and not attached):
[1] Rcpp_0.12.11 compiler_3.4.1 RColorBrewer_1.1-2 plyr_1.8.4
[5] bindr_0.1 tools_3.4.1 digest_0.6.12 tibble_1.3.3
[9] gtable_0.2.0 lattice_0.20-35 pkgconfig_2.0.1 rlang_0.1.1
[13] GGally_1.3.1 parallel_3.4.1 mvtnorm_1.0-6 loo_1.1.0
[17] bindrcpp_0.2 gridExtra_2.2.1 coda_0.19-1 stats4_3.4.1
[21] grid_3.4.1 reshape_0.8.6 inline_0.3.14 glue_1.1.1
[25] R6_2.2.2 rethinking_1.59 magrittr_1.5 scales_0.4.1
[29] codetools_0.2-15 matrixStats_0.52.2 MASS_7.3-47 assertthat_0.2.0
[33] colorspace_1.3-2 labeling_0.3 lazyeval_0.2.0 munsell_0.4.3

reginfotest.csv (7.6 KB)

Bob_Carpenter · August 14, 2017, 3:42pm

There’s a missing semicolon after the 2. By “index position” it means which index.

Look up the line number and see what’s going wrong.

ytan · August 15, 2017, 3:52am

Dear Mr. Carpenter,

thank you for your advice.
A semicolon is located in the next line.

But, after reading your comment, I noticed that an index of the same line was wrong, mu_H_raw[p,t], and changed it into mu_H_raw[p,t-2].

Then, I ran this program, and got the different error message.

SAMPLING FOR MODEL ‘model-treghapSNC_simpmatx2’ NOW (CHAIN 1).
Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: normal_lpdf: Location parameter is nan, but must be finite! (in ‘model39dc68f93d8_model_treghapSNC_simpmatx2’ at line 25)
…
Initialization between (-2, 2) failed after 100 attempts.
Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] “Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.”
[1] “error occurred during calling the sampler; sampling not done”

Stan code

data {
  int P;
  int T;
  int T_pred;
  real<lower=0, upper=10> H[P,T];
}
 
parameters {
  real<lower=0> s_mu;
  real mu_H_raw[P,T-2];
  real<lower=0> s_H;
}
 
transformed parameters {
  real mu_H[P,T];
  for (t in 3:T)
    for (p in 1:P)
      mu_H[p,t] = 2 * mu_H[p,t - 1] - mu_H[p,t - 2] + s_mu * mu_H_raw[p,t - 2]; //reparameterization
 } 
 
model {
  for (t in 3:T)
    for (p in 1:P){
      mu_H_raw[p,t] ~ normal(0, 1); //reparameterization
      H[p,t] ~ normal(mu_H[p,t], s_H);
    }
}
 
generated quantities {
  real mu_all[P,T + T_pred];
  real h_pred[P,T_pred];
  for (t in 1:T)
    for (p in 1:P)
      mu_all[p,t] = mu_H[p,t];
  for (t in 1:T_pred)
    for (p in 1:P){
      mu_all[p,T + t] = 2 * mu_all[p,t - 1] - mu_all[p,t - 2] + s_mu * normal_rng(0, 1);
      h_pred[p,t] = normal_rng(mu_all[p,T + t], s_H);
    }
}

I looked up the line 25, and supposed that mu_H[p,t] had a problem about the value constraint
because H[P,T] can take only integers 0 to 10.

Following my Stan textbook (Matsuura. K, 2016, p168 [in Japanese]), I added two lines in the model block to avoid including unexpected values.

  if (mu_H[p,t] > 10 || mu_H[p,t] < 0 )
  target += negative_infinity();

Stan code

data {
  int P;
  int T;
  int T_pred;
  real<lower=0, upper=10> H[P,T];
}
 
parameters {
  real<lower=0> s_mu;
  real mu_H_raw[P,T-2];
  real<lower=0> s_H;
}
 
transformed parameters {
  real mu_H[P,T];
  for (t in 3:T)
    for (p in 1:P)
      mu_H[p,t] = 2 * mu_H[p,t - 1] - mu_H[p,t - 2] + s_mu * mu_H_raw[p,t - 2]; //reparameterization
 } 
 
model {
  for (t in 3:T)
    for (p in 1:P){
      mu_H_raw[p,t] ~ normal(0, 1); //reparameterization
      if (mu_H[p,t] > 10 || mu_H[p,t] < 0 )
      target += negative_infinity(); // to avoid including unexpected values
      H[p,t] ~ normal(mu_H[p,t], s_H);
    }
}
 
generated quantities {
  real mu_all[P,T + T_pred];
  real h_pred[P,T_pred];
  for (t in 1:T)
    for (p in 1:P)
      mu_all[p,t] = mu_H[p,t];
  for (t in 1:T_pred)
    for (p in 1:P){
      mu_all[p,T + t] = 2 * mu_all[p,t - 1] - mu_all[p,t - 2] + s_mu * normal_rng(0, 1);
      h_pred[p,t] = normal_rng(mu_all[p,T + t], s_H);
    }
}

But, the same error message appeared.

Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: normal_lpdf: Location parameter is nan, but must be finite! (in ‘model39dc38dd7471_model_treghapSNC_simpmatx2’ at line 27)
…
Initialization between (-2, 2) failed after 100 attempts.
Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] “Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.”
[1] “error occurred during calling the sampler; sampling not done”

I also set initial values in the R code as follows after reading this message, but it did not work.

P = nrow(d)
data <- list(P = nrow(d),T = ncol(d) - 1,T_pred = 3,H = d[,-1])
stanmodel <- stan_model(file='model/ch3/model-treghapSNC_simpmatx2.stan')
fit <- sampling(
  stanmodel,
  data = data,
  seed = 243,
  chains=4,
  iter=2000,
  warmup=1000,
  thin=1,
  init=function(){
    list(mu_h[p,1]=rnorm(P,6.2,0.1),s_H=1)
  }
)

dlakelan · August 15, 2017, 4:37pm

It looks to me like you need to initialize mu_h[p,1] mu_h[p,2] to some values before entering this loop:

ytan · August 16, 2017, 11:59pm

Dear Mr.Lakeland,

thank you for your comment.
First, I ran the following codes to initialize mu_h[p,1] mu_h[p,2] .

Stan code

data {
  int P;
  int T;
  int T_pred;
  real<lower=0, upper=10> H[P,T];
}
 
parameters {
  real mu_h1;
  real mu_h2;
  real mu_H_raw[P,T-2];
  real<lower=0> s_mu;
  real<lower=0> s_mu_h1;
  real<lower=0> s_mu_h2;
  real<lower=0> s_H;
}
 
transformed parameters {
  real mu_H[P,T];
  for (t in 3:T)
    for (p in 1:P)
      mu_H[p,t] = 2 * mu_H[p,t - 1] - mu_H[p,t - 2] + s_mu * mu_H_raw[p,t - 2]; //reparameterization
 } 
 
model {
  for (p in 1:P){
  mu_H[p,1] ~ normal(mu_h1,s_mu_h1); // assuming mu_H[p,1] takes similar values across p.
  mu_H[p,2] ~ normal(mu_h2,s_mu_h2); // assuming mu_H[p,2] takes similar values across p.
  }
  for (t in 3:T)
    for (p in 1:P){
      mu_H_raw[p,t] ~ normal(0, 1); //reparameterization
      H[p,t] ~ normal(mu_H[p,t], s_H);
    }
}
 
generated quantities {
  real mu_all[P,T + T_pred];
  real h_pred[P,T_pred];
  for (t in 1:T)
    for (p in 1:P)
      mu_all[p,t] = mu_H[p,t];
  for (t in 1:T_pred)
    for (p in 1:P){
      mu_all[p,T + t] = 2 * mu_all[p,t - 1] - mu_all[p,t - 2] + s_mu * normal_rng(0, 1);
      h_pred[p,t] = normal_rng(mu_all[p,T + t], s_H);
    }
}

R code

setwd("~/docdis")
library(rstan)
library(dplyr)
library(tidyr)
#import data
d <- read.csv(file ='input/ch3/reginfotest.csv', header = TRUE,
              stringsAsFactors = FALSE)

#long -> wide
d <- spread(d, key = year, value = mean)
P = nrow(d)
data <- list(P = nrow(d),T = ncol(d) - 1,T_pred = 3,H = d[,-1])
stanmodel <- stan_model(file='model/ch3/model-treghapSNC_simpmatx3.stan')
fit <- sampling(
  stanmodel,
  data = data,
  seed = 243,
  chains=4,
  iter=2000,
  warmup=1000,
  thin=1,
  init=function(){
    list(mu_h1=6.2, mu_h2=6.2, s_mu_h1=1, s_mu_h2=1, s_H=1)
  }
)

I set initial values based on data. But, I got the following warning and error message.

DIAGNOSTIC(S) FROM PARSER:
Warning (non-fatal):
Left-hand side of sampling statement (~) may contain a non-linear transform of a parameter or local variable.
If it does, you need to include a target += statement with the log absolute determinant of the Jacobian of the transform.
Left-hand-side of sampling statement:
mu_H[p, 1] ~ normal(…)
Warning (non-fatal):
Left-hand side of sampling statement (~) may contain a non-linear transform of a parameter or local variable.
If it does, you need to include a target += statement with the log absolute determinant of the Jacobian of the transform.
Left-hand-side of sampling statement:
mu_H[p, 2] ~ normal(…)

SAMPLING FOR MODEL ‘model-treghapSNC_simpmatx3’ NOW (CHAIN 1).
Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: normal_lpdf: Random variable is nan, but must not be nan! (in ‘model46f05584b36_model_treghapSNC_simpmatx3’ at line 27)
…
Initialization between (-2, 2) failed after 100 attempts.
Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] “Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.”
[1] “error occurred during calling the sampler; sampling not done”

Next, I constrained mu and sigma of the normal distribution to generate mu_H[p,1] and mu_H[p,2] by giving fixed values in the model block.
but the same problem happens.

Stan code

data {
  int P;
  int T;
  int T_pred;
  real<lower=0, upper=10> H[P,T];
}
 
parameters {
  real mu_H_raw[P,T-2];
  real<lower=0> s_mu;
  real<lower=0> s_H;
}
 
transformed parameters {
  real mu_H[P,T];
  for (t in 3:T)
    for (p in 1:P)
      mu_H[p,t] = 2 * mu_H[p,t - 1] - mu_H[p,t - 2] + s_mu * mu_H_raw[p,t - 2]; //reparameterization
 } 
 
model {
  for (p in 1:P){
  mu_H[p,1] ~ normal(6.2,1); // giving fixed values for mu and sigma
  mu_H[p,2] ~ normal(6.2,1); // giving fixed values for mu and sigma
  }
  for (t in 3:T)
    for (p in 1:P){
      mu_H_raw[p,t] ~ normal(0, 1); //reparameterization
      H[p,t] ~ normal(mu_H[p,t], s_H);
    }
}
 
generated quantities {
  real mu_all[P,T + T_pred];
  real h_pred[P,T_pred];
  for (t in 1:T)
    for (p in 1:P)
      mu_all[p,t] = mu_H[p,t];
  for (t in 1:T_pred)
    for (p in 1:P){
      mu_all[p,T + t] = 2 * mu_all[p,t - 1] - mu_all[p,t - 2] + s_mu * normal_rng(0, 1);
      h_pred[p,t] = normal_rng(mu_all[p,T + t], s_H);
    }
}

R code

setwd("~/docdis")
library(rstan)
library(dplyr)
library(tidyr)
#import data
d <- read.csv(file ='input/ch3/reginfotest.csv', header = TRUE,
              stringsAsFactors = FALSE)

#long -> wide
d <- spread(d, key = year, value = mean)
P = nrow(d)
data <- list(P = nrow(d),T = ncol(d) - 1,T_pred = 3,H = d[,-1])
stanmodel <- stan_model(file='model/ch3/model-treghapSNC_simpmatx4.stan')
fit <- sampling(
  stanmodel,
  data = data,
  seed = 243,
  chains=4,
  iter=2000,
  warmup=1000,
  thin=1,
  init=function(){
    list(s_H=1)
  }
)

DIAGNOSTIC(S) FROM PARSER:
Warning (non-fatal):
Left-hand side of sampling statement (~) may contain a non-linear transform of a parameter or local variable.
If it does, you need to include a target += statement with the log absolute determinant of the Jacobian of the transform.
Left-hand-side of sampling statement:
mu_H[p, 1] ~ normal(…)
Warning (non-fatal):
Left-hand side of sampling statement (~) may contain a non-linear transform of a parameter or local variable.
If it does, you need to include a target += statement with the log absolute determinant of the Jacobian of the transform.
Left-hand-side of sampling statement:
mu_H[p, 2] ~ normal(…)

SAMPLING FOR MODEL ‘model-treghapSNC_simpmatx4’ NOW (CHAIN 1).
Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: normal_lpdf: Random variable is nan, but must not be nan! (in ‘model46f030c02b2_model_treghapSNC_simpmatx4’ at line 23)
…
Initialization between (-2, 2) failed after 100 attempts.
Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] “Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.”
[1] “error occurred during calling the sampler; sampling not done”

Bob_Carpenter · August 17, 2017, 1:19am

Was there no other error message for the inits? Have you updated to Stan 2.16? Some of the warnings were getting lost.

The problem here is that random inits are not returning finite log density values. you need to trace which variable is being used before it is defined, or which index it out of bound, etc. That’s why you need the warning messages.

ytan · August 17, 2017, 1:50am

My Rstan version is 2.16.2.
Below the warning message I quoted, I found these lines. In the final part, I found a warning message. Do you mean this?

In file included from C:/Users/yohei/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
from C:/Users/yohei/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
from C:/Users/yohei/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
from C:/Users/yohei/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
from C:/Users/yohei/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
from C:/Users/yohei/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
from C:/Users/yohei/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
from C:/Users/yohei/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
from file249c62e84a86.cpp:8:
C:/Users/yohei/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: “BOOST_NO_CXX11_RVALUE_REFERENCES” redefined

define BOOST_NO_CXX11_RVALUE_REFERENCES

^
:0:0: note: this is the location of the previous definition

Topic		Replies	Views
Divergent transitions in state space modelling with Poisson data Modeling	2	933	June 1, 2017
Latent variables, divergent transitions and prior selection Modeling rstan , techniques , fitting-issues	22	1417	September 19, 2021
Two questions: ①Rejecting initial value but still sampling. ②regarding divergent transitions Modeling fitting-issues	11	1089	November 8, 2020
Help on a state space model Modeling fitting-issues	6	3768	September 3, 2018
Divergent transitions with hierarchical model Modeling	5	722	July 15, 2019

Divergent transitions & BFMI low in a state-space model

Stan code

Stan code

Stan code

R code

Stan code

R code

define BOOST_NO_CXX11_RVALUE_REFERENCES

Related topics