piecewise (0,inf) -> (-inf,inf) transform?

[Bob Carpenter]

As is, when a variable is declared in Stan with
range (0,inf), a log transform is applied to unconstrain
it. In practice, this means we inverse transform it using
the exponential function. This is costly computationally.

Instead, we could use a transform f : (0,inf) -> (-inf,inf)
for some lower bound epsilon > 0 with

f(x) = log x + (1 - log epsilon) if x < epsilon, and

= x / epsilon if x >= epsilon.

Then the inverse function g : (-inf,inf) -> (0,inf) would be

g(y) = exp(y) - (espilon - log epsilon) if y < log epsilon, and

= epsilon * y if y >= log epsilon.

The function g(y) is continuous and differentiable everywhere with

g'(y) = exp(y) if y < log epsilon

= epsilon if y >= log epsilon

leading to log Jacobian

log | g'(y) | = y if y < log epsilon, and

= log epsilon if y >= log epsilon.

This will lead to faster transforms and fewer derivative terms
on the auto-dif stack. It also prevents overflow as y -> inf.

I'm not clear on whether this will be a good or bad thing
geometry-wise for HMC.

- Bob

[Michael Betancourt]

I've played a fair number of games with these kinds of transforms and I've never found a satisfactory alternative to exp/log transforms. The problem in your transform, for example, is that, while the transform and its first derivative are continuous, its second derivative (and hence in some sense the optimal step size) is not. I'll have to take the opportunity to advocate for bouncing again -- it removes the need for the expensive transform and does not introduce any new Jacobians.

-Mike