Unexplained poor sampling for simple state space model involving a cyclic component

I do not succeed to get reasonable sampling when a cyclic component is added to a state-space model. The same model and data using PyMC provides excellent sampling. Maybe that my code for a cyclic component is wrong/sub-optimal. Any suggestion would be appreciated.
The code below implements a combination of integrated random walk with a cycle of known frequency.

data{
  int<lower=1>    N;    // Time series length
  vector[N]       y;       // Observations
  real            omega; // normalized angular frequency [rad] of the cycle
}

parameters{
  real slope_ini;
  real level_ini;
  real cycle_ini;
  real conj_cycle_ini;
  vector[N-1]  slope;
  vector[N]  cycle;
  vector[N]  conj_cycle;
  real<lower=0>   sigma_slope; 
  real<lower=0>	  sigma_irreg;
  real<lower=0>   sigma_cycle;
}

transformed parameters {
  vector[N] level;
  level[1] = level_ini; 
  for (n in 2:N) {
    level[n] = level[n-1] + slope[n-1];
  }  
}

model{
  // observation equations
  y[1:N] ~ normal(level[1:N] + cycle[1:N], sigma_irreg);

  // slope equations
  slope[1] ~ normal(slope_ini, sigma_slope);
  for (n in 2:N-1) {
    slope[n] ~ normal(slope[n-1], sigma_slope); 
  }  
    
  // cycle equations
  cycle[1] ~ normal(cycle_ini, sigma_cycle);
  conj_cycle[1] ~ normal(conj_cycle_ini, sigma_cycle);  
  for (n in 2:N) {
    cycle[n] ~ normal(cos(omega)*cycle[n-1] + sin(omega)*conj_cycle[n-1], sigma_cycle);
    conj_cycle[n] ~ normal(-sin(omega)*cycle[n-1] + cos(omega)*conj_cycle[n-1], sigma_cycle);
  } 	
}

```stan

Can you be more specific? Were there divergence warnings or low effective sample size or multiple chains failing to converge?

It would help if you could share the working PyMC code and the calls used to fit it. Were you using anything in PyMC or Stan other than random initialized NUTS? If you were just using PyMC’s NUTS implementation, my guess would be that it’s not the same model being implemented.

The main issue I see here computationally is that you’re using a centered parameterization of a hierarchical model for cycle and conj_cycle. My guess is that your program will mix much better with a non-centered parameterization, which only requires a change in declaration (and a recent enough version of Stan).

vector<offset= cycle_ini, multiplier=sigma_cycle>[N] cycle;

You can speed up that Stan model computation. Here are some suggestions.

  1. Use cumulative sum to define level in terms of slope.
level[1] = level_ini;
level[2:] = level_ini + cumulative_sum(slope);

or as a one-liner with a compound declare/define.

vector[N] level = level_ini + append_row(0, cumulative_sum(slope));
  1. Full indexes can be dropped, so what you’re calling “observation equations” can be just:
y ~ normal(level + cycle, sigma_irreg);

This will be more efficient because it doesn’t do an allocate and copy for y, level, and cycle.

  1. The loop in the model block can be replaced with vectorized forms:
cycle[2:] ~ normal(cos(omega) * cycle[1:(N - 1)] + sin(omega) * conj_cycle[1:(N-1])]);

and same for conj_cycle. It’d be more efficient still to define cycle and conj_cycle to be of length N-1 so this doesn’t require any slicing and then cycle[1] and conj_cycle[1] can be scalar parameters.

  1. Precompute cos(omega) and sin(omega) and store as transformed data variables and reuse.
transformed data {
  real cos_omega = cos(omega0;
  real sin_omega = sin(omega);
  real neg_sin_omega = -sin_omega;
}

Dear Bob,

All convergence measures failed for the slope disturbance variance: poor sampling, strong autocorrelations over hundreds of lags, R-hat etc. Mixing very poor. I also found that even a simple integrated random walk model shows similar problems with the slope disturbance variance. Only the observation disturbance variance converges. Using the same (synthetic) dataset and model in PyMC gives perfect results. Therefore, I suspect that my model implementation is at least sub-optimal. I will first follow your recommendations regarding the simple integrated random walk model and will report the results. If I get this simple model working I will add a cycle, follow your instructions again, and let you know.

Anyway thanks a lot for your response.

Regards, Roland Klees

I converted the jupyter notebooks to Python. Enclosed the two files.

Roland Klees

SSM_cmdstanpy-experimental_synthetic_dataset.py (13.4 KB)

SSM_pymc-experimental_synthetic_dataset.py (26.8 KB)