bob-carpenter
commented
2 years ago To answer you question, if the transforms and adjoint-Jacobian product are expressed as functions, then you can pull their generated C++ out and do that test in C++ using Stan (there's a Jacobian functional, absolute value function, and log determinant function). I would instead suggest adding a few test cases where we know the answer and instead test that the log density computed by Stan is correct. That we can do without opening C++. We can also test that we get the right answer for sampling, which won't be true if the Jacobian adjustment is wrong. Of course, the Jacobian might be subtly wrong, and we probably wouldn't catch that with a probabilistic test like MCMC.
So I guess what I'm saying is that I'd want to do those tests exactly as @sethaxen proposed for C++ production code, but I think it's overkill for this paper. Having said that, I can help show how to do it if you're still interested (the main pain point is getting everything made, so it'd be easier to add it as a standard Stan unit test, so the dev ops is all in place).
I've been locally testing my Jacobian derivations in Julia with AD, but ideally our Stan implementations would be tested for correctness as well by computing the Jacobian of the actual transforms using AD and comparing its log-abs-det with the log log-abs-det-jac terms returned by the implementations. Is this straightforward to do? This is outside of my Stan expertise.