When sampling bounded variables using NUTS, does pyMC first sample in terms of a continuous variable transform, allowing for well behaved calculations of the gradient, and then map these samples to the bounded domain?
If so when using my own custom likelihood function with bound random variable priors, should the gradient of my custom likelihood be in terms of the continuous variable transform or the untransformed variable?
Yes. You just have to make sure the appropriate default transform is being used, by either passing it manually when you create a model variable or making your distribution inherit from classes like PositiveContinuous, UnitContinuous, which specify a default transformation.
Should be in terms of the bounded/constrained space. The adjustment when going from unbounded to bounded space is done by PyMC and is transform specific.