# Numerical methods for calculating derivative of gamma function

I am working a mathematical statistic problem. It requires me to find MLE for a Gamma distribution by using numerical methods (rather than Method of Moments). Because Gamma distribution has two parameters, I have to solve a system of two non-linear equations which involve gamma function and its derivatives. I know there are methods that approximate gamma functions, Spouge's approximation or Lanczos approximation. (This post has a good summary Algorithm to compute Gamma function.) But I don't know how to approximate the derivative of a gamma function, that is, $$\Gamma'(z) = \int_0^\infty (z-1) t^{z-1} e^{-t} \ln t \; dt .$$ I guess one can use Riemann sum to approximate this value, since the tail when $t \to \infty$ must be bounded. But this may be inaccurate. Is there any better way to approximate this integral (derivative of gamma function). Thanks a lot!

• In Matlab you could do this with the functions psi and gamma: $\Gamma'(x)=\psi(x) \Gamma(x)$. – Ian Feb 4 '17 at 15:20
• Hi, Ian. Thank you! I haven't used MATLAB before. I will check it definitely. I have one question. How good are gamma and psi performance! – jwyao Feb 4 '17 at 15:45
• Why not just try the benchmarking yourself? It's not like the programming is difficult. – Ian Feb 4 '17 at 15:47
• Sure. I am reading MathWorks document. It seems MATLAB can do it. Thank you so much! – jwyao Feb 4 '17 at 15:49