_ln_dgamma(x, a, b) = a*log(b) - lgamma(a) + (a-1)*log(x) - b*x
pgamma(x, shape, rate) = (x<0)? 0 : igamma(shape, x*rate)
dgamma(x, shape, rate) = (x<0)? 0 : (x==0)? ((shape<1)? 1/0 : (shape==1)? rate : 0) : (rate==0)? 0 : exp(_ln_dgamma(x, shape, rate))
gamma_offset = 0.0
gamma_k = 15.7880378786862
gamma_theta = 0.980078674752871
gamma_fit(x) = 100.0 * dgamma(x - gamma_offset, gamma_k, gamma_theta)
#gamma_dist(x) = x**(gamma_k - 1) * exp((-1.0) * x / gamma_theta) / (gamma_theta**gamma_k * gamma(gamma_k))
#gamma_fit(x) = 100.0 * gamma_dist(x - gamma_offset)