Unexpected floats generated

With these inputs

{"J":2,"N":4,"D":2,"K":3}

this program

data {
  int<lower=0> N;
  int<lower=1> J;
  int<lower=1> K;
}
generated quantities {
  real x[N];
  real u[N];
  for (n in 1:N) {
    x[n] = normal_rng(0,1);
    u[n] = normal_rng(0,1);
  }
  real beta[K];
  real gamma[K];
  vector[K] b0;
  vector[K] g0;
  for (k in 1:K) {
      beta[k] = normal_rng(1,1);
      gamma[k] = normal_rng(3,1);
      b0[k] = normal_rng(0,1);
      g0[k] = normal_rng(0,1);
  }
  real<lower=0> sigma = uniform_rng(.6,.8);
  real mu[N];
  real y[N];
  simplex[K] theta;
  int<lower=1,upper=J> jj[N];
  int<lower=1,upper=K> kk[N];
  theta = dirichlet_rng(rep_vector(1,J));
  vector[K] h0[J];
  for (j in 1:J) {
      for (k in 1:K) {
        h0[j][k] = normal_rng(0,1);
      }
  }
  for (n in 1:N) {
    jj[n] = categorical_rng(theta);
    kk[n] = categorical_logit_rng(h0[jj[n]] + x[n]*b0 + u[n]*g0);
    mu[n] = x[n] * beta[kk[n]] + u[n] * gamma[kk[n]];
    y[n] = normal_rng(mu[n], sigma);
  }
}

produces a csv whose first two lines are

lp__,accept_stat__,x.1,x.2,x.3,x.4,u.1,u.2,u.3,u.4,beta.1,beta.2,beta.3,gamma.1,gamma.2,gamma.3,b0.1,b0.2,b0.3,g0.1,g0.2,g0.3,sigma,mu.1,mu.2,mu.3,mu.4,y.1,y.2,y.3,y.4,theta.1,theta.2,theta.3,jj.1,jj.2,jj.3,jj.4,kk.1,kk.2,kk.3,kk.4,h0.1.1,h0.2.1,h0.1.2,h0.2.2,h0.1.3,h0.2.3
0,0,-0.980565,1.64421,-0.0895129,0.448958,0.143304,-0.971958,-0.273569,-0.396082,0.667543,1.15117,1.30875,2.3475,4.07715,1.65089,0.420343,0.967807,0.954116,-1.56345,-0.116646,0.930609,0.769948,-1.04674,-1.1841,-0.701959,-0.630105,-1.56455,-1.97479,-0.24165,-1.48083,0.579633,0.420367,1,2,2,1,3,1,1,1,1.66301,-0.424079,0.0171969,-0.0980869,1.59318,0.101672,0

containing as kk.4 floating-point values where I expected small integers. How can I fix it?

Looking at the csv, it looks like all of the outputs are off by 1 (i.e., theta.3 appears to have the integer value for jj.1. Did you have initialisation errors about NaN or theta not being a simplex when you ran this?

I believe the cause is your construction of the rng for theta:

  simplex[K] theta;
...
  theta = dirichlet_rng(rep_vector(1,J));

You’ve declared theta as a simplex of size K, but your dirichlet rng is returning a vector of size J. Given the data input that you posted above, this means that the first two elements of theta are initialised to a two-element simplex, and the last element is uninitialised (i.e., NaN)

No, no message like that, but that did solve the issue.

1 Like