Vectorization of the normal_rng and bernouli_rng

Modeling
1

Hi there, my generated quantities block needs to be vectorized and I appreciate if anyone could help me with that. Here is the code that needs to be vectorized:

    generated quantities {      
      vector[N_pred] y_pred;
      for (n in 1:N_pred){
        if (bernoulli_logit_rng(X[n] * beta1)) { //beta1 is defined as vector[I], X is a matrix of N_pred * D
          y_pred[n] = normal_rng(X[n] * beta2, sigma); //beta2 is defined as vector[I], sigma is a real value
       }
      }
   }

The following attempts give me error as below:

    generated quantities {      
      vector[N_pred] y_pred= normal_rng(X * beta2, sigma);
      vector[N_pred] labels = bernoulli_logit_rng(X * beta1);
      y_pred = y_pred  .* labels;
   }

and the error is :

ValueError: SYNTAX ERROR, MESSAGE(S) FROM PARSER:

No matches for: 

  normal_rng(vector, real)

Available argument signatures for normal_rng:

  normal_rng(real, real)

I tried changing sigma to rep_vector(1, N_pred) * sigma to have it as normal_rng(vector, vector) but still get an error. I believe I am doing similar thing to this, so I can not figure out what is wrong. I appreciate if anyone could take a look please.

PS: I am using CmdStan.

2

I’ve used this approach:
vector[K] beta = to_vector(normal_rng(rep_vector(0,K), rep_vector(1,K)));

3

@maswiebe but it still gives “No matches for normal_rng(vector, vector)”. right? Previously I tried rep_vector(1, N_pred) * sigma to have both mean and var of normal as vectors but it did not work.

4

The issue here is with the return type actually. The vectorised rng functions will always return an array, rather than a vector. So you just to convert the result to a vector before you assign it to an output:

generated quantities {      
      vector[N_pred] y_pred = to_vector(normal_rng(X * beta2, sigma));
      vector[N_pred] labels = to_vector(bernoulli_logit_rng(X * beta1));
      y_pred = y_pred  .* labels;
}
4 Likes
5

Very helpful, @andrjohns. I was just pulling my hair out over this, and your comment saved the day.

2 Likes
6

Hey @andrjohns, is there a way to build on these sensibilities to compute the log_lik with vectorized style syntax, instead of the typical for loop syntax like shown in the loo documentation?

From the loo documentation, we find:

generated quantities {
  vector[N] log_lik;
  for (n in 1:N) {
    log_lik[n] = bernoulli_logit_lpmf(y[n] | X[n] * beta);
  }
}
7

Afraid not. We’ve been meaning to make a version of our densities that returns an array but we haven’t yet.

The way to make this a bit faster is to lean into the matrix arithmetic as follows:

vector[N] logit_p = X * beta;
for (n in 1:N) { 
  log_lik[n] = bernoulli_logit_lpmf(y[n] | logit_p[n]);
}

The matrix arithmetic is a lot faster than pulling out row vectors and multiplying.

I wish that was a problem I could have.

3 Likes
8

Got it. Thanks for the reply, @Bob_Carpenter.