random log-gamma for taming underflow issues

[chelate]

Hi all sorry about this. I just have this function from an old project, it works really well for log-gamma sampling for small shape parameter

function randlogGamma(a;scale = 1)
    if a > .2
        return log.(rand(Gamma(a))) #fall back is fine, just not perfect Float-dense
    else
        L = 1 / a - 1
        w = a / MathConstants.e / (1 - a)
        ww = 1 / (1 + w)
        η = z -> ( z >= 0 ? exp(-z) : w * L * exp(L * z))
        h = z -> exp(-z - exp(-z / a))
        function rh(a)
            while true
                U = rand()
                z = (U <= ww) ? -log(U / ww) : log(rand()) / L
                h(z) / η(z) > rand() && return(z)
            end
        end
        return log(scale) - rh(a) / a
    end
end

It is especailly useful for stable beta sampling when α/β are super small. I don't know who wrote the code anymore, but it's very useful, I've had to copy it to a few projects now.

The problem is, without this you don't get the right power laws because of underflow for small a. Ideally, for best possible numerical stability, beta and Dirichlet would always be defined through logsumexp acting on vectors of these guys to enforce normalization and keep log(0) at bay.

I know just posting code instead of a pr is terrible manners and I will see if I can figure out how to do that. But we Really need a log-gamma sampler to do gamma and beta properly!

[chelate]

see #1003 #1024 also found the paper where I got it from: https://arxiv.org/pdf/1302.1884.pdf