Passing matrix with one dimension equal to 0 works only sometimes?

Modeling
1

Why can I pass a matrix with one dimension equal to 0 to Stan in one case but not in another? The example below works (but shouldn’t?) while in another case it unapologetically states that

  Exception: mismatch in number dimensions declared and found in context; processing stage=data initialization; variable name=x_e; dims declared=(150,0); dims found=(0)  (in 'model_msm_constant_continuous_v2' at line 15)

Working example:

library(rstan)

m <- rstan::stan_model("test.stan")

N <- 300
K <- 0
x <- matrix(rbinom(N * K, 1, 0.3), ncol = K, nrow = N)
alpha <- rnorm(1, 0, 1)
beta <- rnorm(K, 0, 1)
sigma <- abs(rnorm(1, 0, 1))
y <- rnorm(N, alpha +  x %*% beta, sigma)


data <- list(
  N = N,
  K = K,
  x = x,
  y = y
)

fit <- rstan::sampling(m, data = data, chains = 2, iter = 1000)

data {
  int N;
  int<lower = 0> K;
  matrix[N, K] x;
  vector[N] y;
}

transformed data {
  int K_ = K == 0 ? 1 : K;
  matrix[N, K_] x_ = K == 0 ? rep_matrix(0.0, N, K_) : x;
  print(K);
  print(x);
}

parameters {
  real alpha;
  vector[K_] beta;
  real<lower = 0> sigma;
}

model {
  target += normal_lpdf(beta | 0, 1);
  target += normal_lpdf(y | alpha + x_ * beta, sigma);
}

This breaks if Mx_e is equal to 0 and x_e[150, Mx_e]

data {
  // dimensions and slicing
  int<lower = 1> N;                 // observations
  int<lower = 2> T;                 // time points
  int<lower = N * T> NT;            // product
  int startstop[N, 2];              // slicer
  int<lower = 2> K;                 // N(states)

  // predictor dimensions
  int<lower = 0> Mx_d;            // N(varying parameters of continuous process)
  int<lower = 0> Mx_e;            // N(varying parameters of continuous process)

  // predictors
  matrix[NT, Mx_d] x_d;         // matrix of shared continuous of discrete process
  matrix[NT, Mx_e] x_e;         // matrix of varying continuous of discrete process

  // output
  vector[NT] y;                     // vector of output

  // intercepts
  int has_intercept[5];             // 1: alpha, 2: beta

  // shared
  int shared_TP;                       // 0: individual, 1: shared
}

transformed data {
  // add intercepts
  int Mx_d_ = has_intercept[1] ? Mx_d + has_intercept[1] : (Mx_d == 0 ? 1 : Mx_d);
  int Mx_e_ = has_intercept[2] ? Mx_e + has_intercept[2] : (Mx_e == 0 ? 1 : Mx_e);
  matrix[NT, Mx_d_] x_d_ = has_intercept[1] ? append_col(rep_vector(1.0, NT), x_d):
                                             (Mx_d == 0       ? rep_matrix(0.0, NT, 1) : x_d);
  matrix[NT, Mx_e_] x_e_ = has_intercept[2] ? append_col(rep_vector(1.0, NT), x_e):
                                             (Mx_e == 0       ? rep_matrix(0.0, NT, 1) : x_e);

  int N_ = shared_TP ? 1 : N;
}
1 Like
2

‘solved’ A matrix with 0 columns and >0 rows is apparently treated as a matrix of NAs and thereby empty. Creating that empty matrix with any value fixes the problem, e.g.

x_d <- matrix(0, ncol = 0, nrow = NT)

works but

x_d <- matrix(NA, ncol = 0, now = NT)

does not despite them looking the same when called.

3 Likes
3

Fascinating!

> a_na <- matrix(NA, nrow = 0, ncol = 10)
> a_na
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
> str(a_na)
 logi[0 , 1:10] 
> b_na <- matrix(0, nrow = 0, ncol = 10)
> b_na
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
> str(b_na)
 num[0 , 1:10]

So one of them is a logical, whilst the other a numeric.

2 Likes
4

Good observations! Just to note that R has NA_real_ for when you want to specify NAs that are numeric:

a_na <- matrix(NA_real_ , nrow = 0, ncol = 10)
str(a_na)
# num[0 , 1:10]
1 Like