How to vectorize a for loop in Stan?

Modeling
1

Apologies for the simple question, but I’m having trouble understanding how vectorization works in Stan.

The following code with a for loop works correctly

for (n in 1:N) {
  target += log(lognormal_cdf(y_upper[n] | mu, sigma) - lognormal_cdf(y_lower[n] | mu, sigma));
}

while the “vectorized” code I thought was equivalent does not:

target += log(lognormal_cdf(y_upper | mu, sigma) - lognormal_cdf(y_lower | mu, sigma));

Could someone explain how to correctly vectorize the for loop?

Full code

data {
  int<lower=0> N;            
  array[N] int y_lower;      
  array[N] int y_upper;   
}
parameters {
  real mu;
  real<lower=0> sigma;
}
model {
  for (n in 1:N) {
  target += log(lognormal_cdf(y_upper[n] | mu, sigma) - lognormal_cdf(y_lower[n] | mu, sigma));
  }
}
2

Vectorized probability functions in Stan (as a general rule) return the sum of the element wise application. So the two programs you wrote are not really equivalent, since they change the order of operations

The for loop is computing sum(log(difference of two values)). The statement is computing log(sum(all first values) - sum (all second values)).

For this kind of computation there is not really anything better than the for loop I believe

3 Likes
3

Stan has a vectorized erf, so it should be possible to code a vectorized element-wise lognormal_cdf at the Stan level, but I suppose that’s unlikely to be efficient compared to using Stan’s lognormal_cdf inside a for-loop, since Stan’s implementation is presumably optimized and autodiffed at the C++ level.