I am trying to specify a box-cox model and generate predictions in the original scale.

I get several errors with this message: The following variables have undefined values
The problem is in this line: y_pred[n] = (normal_rng(alpha + log(x[n])*beta, sigma)* lambda + 1)^(1/lambda);

Hmm something funky is going on here. The generated quantities is generating NaNs in the output and Rstan is complaining that nothing was assigned (by default variables declared in generated quantities are NaN). To debug this, throw in some print statements in your generated quantities:

There are some cases in which the normal_rng function result in values that after transforming end being NAs.

Using these values (from printing some iterations):

alpha = -1.26999
logx = 7.90211
beta = 0.101912
lambda = 1.12192
sigma = 0.316502
a = rnorm(100, alpha + beta*logx, sigma)
df = data.frame(value = a, tf = (a * lambda + 1)^(1/lambda))
df[is.na(df$tf),]
value tf
9 -1.1070863 NaN
10 -1.1268851 NaN
16 -1.1499958 NaN
35 -1.2567046 NaN
38 -0.9276948 NaN
39 -0.9124289 NaN
42 -0.9351521 NaN
59 -1.1140359 NaN
89 -1.4590619 NaN

I should limit in some way the lower value I get from the random function.

It looks like in those models they have avoided doing the posterior predictives that youâ€™re using. Maybe for the same reason. Does this seem familiar to you @Stephen_Martin?

Thanks! Yes, I understand. I am trying to find a good solution. I plan to combine several models (stacking) where box-cox is one of them, so I would need predictions to weight them and use them later.

Mind you, I had done a â€śx - min(x) + 1â€ť transform prior to using the stan model, hence the extra little â€ś-1 + stan_data_t$minYâ€ť tidbit.

That seems to have worked.
The gist of it is - Fit the Box-cox model; obtain random draws ON THE BC-TRANSFORMED scale; then in R un-transform those yreps per-iteration. It handles those problems a bit better. You may still get some Inf/NaN problems, mind you, but such is life.

That lambda == 0 case shouldnâ€™t arise unless you explicitly initialize at zero because itâ€™s a measure zero event (not quite with floating point, but close enough). Were you running into problems with that otherwise?

Itâ€™s worth doing the power calculation in a loop, then you can vectorize the matrix product and normal; and you can pull the Jacobian out (Iâ€™m assuming thatâ€™s what it isâ€”I always have to work powers out from first principles) into a single vectorized operation:

Iâ€™m sure the lambda == 0 case never actually occurs, it was included to be â€ścorrectâ€ť. Basically, that code was being used for a simulation paper - Correctness was valued over efficiency, b/c I wanted to avoid a reviewer saying â€śthat BC transform is incomplete! What about lambda = 0â€™s case?â€ť

And youâ€™re right; the code could be made more efficient. Again - correctness (and readability) were valued over efficiency in this particular case.

Thanks @Stephen_Martin and @Bob_Carpenter. I was using a much more crude and probably wrong solution: just to skip those predictions with no valid values and assign a value after 5000 iterations. The uncertainty of the estimates, specially at the bottom of the distribution, becomes weird:

i = 1;
while ((is_nan(y_pred[n])) && (i <= 5000)) {
y_pred[n] = ((normal_rng(alpha + log(x[n])*beta, sigma)) * lambda + 1)^(1/lambda);
i = i + 1;
}
if (is_nan(y_pred[n])) {
y_pred[n] = 0; // very crude solution!
}

I have further questions.

Does the way you compute log_lik make possible to compare with models that do not change the dependent variable using, for instance, WAIC or LOO?

I donâ€™t get well your definition of target. Could you give some hints on that?