weird behavior of exponential_log in model vs generated quantities block; something to do with constraints?

[Seth Flaxman]

Hi. I can't figure out what is going on with this very simple example
(Stan inline w/ R code, apologies in advance):

fit = stan(model_code='
parameters {
real<lower=0> y;
}
model {
increment_log_prob(exponential_log(y,1)); // equivalently: y ~ exponential(1);
}
generated quantities {
real lp;
lp <- exponential_log(y,1);
}')

out=extract(fit)
plot(out$lp__,out$lp)

I'm expecting lp and lp__ to be the same, but they're not. lp gives
exactly what I expect, same as dexp(...) in R (since the rate is 1,
the log-likelihood is -y). Can anyone explain what's going on with
lp__?

For comparison, this works as expected:

fit = stan(model_code='
parameters {
real y;
}
model {
increment_log_prob(normal_log(y,1,1));
}
generated quantities {
real lp;
lp <- normal_log(y,1,1);
}')
out=extract(fit)
plot(out$lp__,out$lp)

And most puzzingly, this does what I expect (although a bunch of
exceptions are thrown whenever y is negative):

fit = stan(model_code='
parameters {
real y;
}
model {
increment_log_prob(exponential_log(y,1));
}
generated quantities {
real lp;
lp <- exponential_log(y,1);
}')
out=extract(fit)
plot(out$lp__,out$lp)

So I guess the behavior in the first model has something to do with
the constraints on y...but that doesn't make sense to me.

Thanks,
Seth

[Bob Carpenter]

Yes, you're getting the Jacobian term from the
inverse of the log-transfrom of y to the unconstrained
space. Specifically, that's an extra log(y) term,
because y is transformed to log(y) and the inverse transform
is exp, and the log derivate of that evaluated at log(y)
is just log(y).

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>

[Daniel Lee]

In the manual, see the chapter on "Transformations of Constrained Variables" for gory details.

To fix your example:
fit = stan(model_code='
parameters {
  real<lower = 0> y;

}
model {
  increment_log_prob(normal_log(y,1,1));
}
generated quantities {
  real lp;
  lp <- normal_log(y,1,1);

  lp <- lp + log(y);
}')
out=extract(fit)
plot(out$lp__,out$lp)