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

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

‘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

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

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