Polygamma Functions: The Trigamma and Invdigamma Functions There seems to be very few implementations of these functions in modern languages that have clear licenses. So I'm making them available here under the MIT license.

invdigamma.jl

# Implementation of fixed point algorithm described in
#  "Estimating a Dirichlet distribution" by Thomas P. Minka, 2000
function invdigamma(y::Float64)
    # Closed form initial estimates
    if y >= -2.22
        x_old = exp(y) + 0.5
        x_new = x_old
    else
        x_old = -1.0 / (y - digamma(1.0))
        x_new = x_old
    end

    # Fixed point algorithm
    delta = Inf
    iteration = 0
    while delta > 1e-12 && iteration < 25
        iteration += 1
        x_new = x_old - (digamma(x_old) - y) / trigamma(x_old)
        delta = abs(x_new - x_old)
        x_old = x_new
    end

    return x_new
end

trigamma.jl

# Trigamma function
#
# Implementation of algorithm described in
# "Algorithm AS 121: Trigamma Function" by B. E. Schneider, 1978
function trigamma(x::Float64)
    trigam = 0.0
    z = x
    if x <= 0.0
        throw(DomainError())
    end

    if x <= 1e-4
        return 1.0 / (z * z)
    end

    while z < 5.0
        trigam += 1.0 / (z * z)
        z += 1.0
    end

    y = 1.0 / (z * z)
    return trigam +
           0.5 * y +
           (1.0 +
            y * (1.0 / 6.0 +
            y * (-1.0 / 30.0 +
            y * (1.0 / 42.0 +
            y * -1.0 / 30.0)))) / z
end