Calculation of derivatives by deriv()

[Chi-lin Tsai]

Hello,

I'm learning deriv(). I use the log gamma function in practice, and use Mata's digamma() and trigamma() to check whether I do it right.

As can be seen from the output below, the first derivative from deriv() is the same as that from digamma(), but the second derivative from deriv(D,2) is very different from what trigamma() produces. I wonder whether the difference is due to numerical errors. If it is not, then I must misunderstand something. Therefore, I would like to seek your helps.

cap m: mata drop D
cap m: mata drop DLnGamma()
mata:
    void DLnGamma(x,y){
        y=lngamma(x)
    }
    D = deriv_init()
    deriv_init_evaluator(D, &DLnGamma())
    deriv_init_params(D, 2)
    Deriv =deriv(D, 1),
           deriv(D, 2)
end
mata: Deriv\(digamma(2),trigamma(2))
                 1             2
    +-----------------------------+
  1 |  .4206921333   5.12566e+11  |
  2 |  .4227843351   .6449340668  |
    +-----------------------------+

Below is a bit more test. I use deriv() to calculate the first derivative of the digamma function. If I understand correctly, the result should produces the second derivative of the log gamma function (i.e., the result from trigamma().) This time, the results do match. However, I still can't figure out why the code above did not produce the second derivative of the log gamma function, while the code below could.

cap m: mata drop D
cap m: mata drop DDiGamma()
mata:
    void DDiGamma(x,y){
        y=digamma(x)
    }
    D = deriv_init()
    deriv_init_evaluator(D, &DDiGamma())
    deriv_init_params(D, 2)
    Deriv =deriv(D, 0),
           deriv(D, 1)
end
mata: Deriv\(digamma(2),trigamma(2))
                 1             2
    +-----------------------------+
  1 |  .4227843351   .6449340674  |
  2 |  .4227843351   .6449340668  |
    +-----------------------------+

Thank you very much in advance.