2

In my Stan code, I want to add an ICAR-term (phi) to the following covariate model:

// covariate models with logit link
  vector[total_surveys] logit_p = (X_p * beta_p) + phi[ii_sampled]; 

However the dim(X_p) = (1900, 3) and dim(beta_p) = 3

Thus dim(X_p * beta_p = 1900

While the dim(phi[ii_sampled]) is 95

Therefore, I would like to do:

phi_p <- rep(phi[ii_sampled], each = 20)

In essence, my problem comes down to doing (example with different dimensions):

phi <- c(1.2, -0.5, 2.1, -0.7)
phi_p <- rep(phi, each = 3)

phi_p
>(1.2, 1.2, 1.2, -0.5, -0.5, -0.5, 2.1, 2.1, 2.1, -0.7, -0.7, -0.7)

in Stan

Unfortunately, the rep() function is not available in Stan. You would have to do a loop instead.

What would such a loop look like?

1 Answers1

3

The following is simple and efficient in Stan, because loops in Stan get compiled down to efficient C++ loops. This could be made shorter, but it'd be less efficient and probably less clear.

functions {
  vector rep_each(vector x, int K) {
    int N = rows(x);
    vector[N * K] y;
    int pos = 1;
    for (n in 1:N) {
      for (k in 1:K) {
        y[pos] = x[n];
        pos += 1;
      }
    }
    return y;
  }
}
transformed data {
  print("rep_each([1, 10, 100]', 4) = ",
        rep_each([1, 10, 100]', 4));
}

Running this prints

rep_each([1, 10, 100]', 4) = [1,1,1,1,10,10,10,10,100,100,100,100]
Bob Carpenter
  • 3,613
  • 1
  • 20
  • 13