/***************  Loggamma Derivatives *****************************
Author:
Paul L. Fackler     Department of Agricultural and Resource Economics
North Carolina State University - Box 8109     Raleigh, NC, USA 27695
919-515-4535                                    paul_fackler@ncsu.edu
Date: May 1998

Code provided gratis without performance guarantees for private
non-commercial use.
*********************************************************************/

@ Calculates the value of the psi (digamma) function, the first
  derivative of the loggamma function.  Accepts positive matrices.
  From Abromowitz and Stegun. Uses formula 6.3.18 with recurrence
  formula 6.3.5. Requires workspace at least 5 times the size of X.
@
proc psi(x);
  local p;
  x=x+6;
  p=1/(x.*x);
  p=(((0.004166666666667*p-0.003968253986254).*p+
      0.008333333333333).*p-0.083333333333333).*p;
  p=p+ln(x)-0.5/x-1/(x-1)-1/(x-2)-1/(x-3)-1/(x-4)-1/(x-5)-1/(x-6);
retp(p); endp;


@
  3 procs to calculate the value of the trigamma, tetragamma, and
  pentagamma functions, the second, third, and fourth derivatives of
  the loggamma function. Accepts positive matrices. From Abromowitz
  and Stegun. Uses formulas 6.4.11 and 6.4.12 with recurrence formula
  6.4.6. Each requires workspace at least 5 times the size of X.
@
proc trigamma(x);
  local p,i;
  x=x+6;
  p=1/(x.*x);
  p=(((((0.075757575757576*p-0.033333333333333).*p+0.0238095238095238)
     .*p-0.033333333333333).*p+0.166666666666667).*p+1)./x+0.5*p;
  i=0; do until i>5;i=i+1;
    x=x-1;
    p=1/(x.*x)+p;
  endo;
  retp(p);
endp;

proc tetragam(x);
  local p,i;
  x=x+6;
  p=1/(x.*x);
  p=(((((-0.833333333333333*p+0.3).*p-0.1666666666666667).*p
     +0.1666666666666667).*p-0.5).*p-1).*p;
  p=-1/(x^3)+p;
  i=0; do until i>5;i=i+1;
    x=x-1;
    p=-2/(x^3)+p;
  endo;
  retp(p);
endp;

proc pentagam(x);
  local p,i;
  x=x+6;
  p=1/(x.*x);
  p=(((((10*p-3).*p+1.3333333333333333).*p-1).*p+2).*p+2+3/x).*p./x;
  i=0; do until i>5;i=i+1;
    x=x-1;
    p=6/(x^4)+p;
  endo;
  retp(p);
endp;

@
  Calculates the value of the polygamma function, the Nth derivative
  of the loggamma function. Accepts positive matrices for X and
  integer matrices for N. X and N must be ExE conformable. From
  Abromowitz and Stegun. Uses formulas 6.4.11 and 6.4.12 with
  recurrence formula 6.4.6.
@
proc polygam(x,n);
  local p,q,i,b;
  x=x+6;
  p=1/(x.*x);
  q=0;
  let b=0 0.166666666666667 0 -0.033333333333333 0 0.023809523809524 0
-0.03333333333333 0 0.075757575757576 0 -0.253113553113553 0
1.16666666666667;
  i=14; do until i<1.5;
    q=(q+b[i,1]*((i+n-1)!/i!)).*p;
    i=i-2;
  endo;
  p=(-1)^(n-1);
  q=p.*((n-1)! + 0.5*(n!)./x + q)./x^n;
  p=p.*n!;
  i=0; do until i>5;i=i+1;
    x=x-1;
    q=q+p./(x^(n+1));
  endo;
  retp(q);
endp;

@
  Calculates the inverse of the PSI function using the method of
  bisection. Uses the fact that the exp(p) < psiinv(p) < exp(p)+1.
@
proc psiinv(p);
  local l,x;
  l=1;
  x=exp(p);
  do until l<1e-8;
    x=x+l*sign(p-psi(x));
    l=l/2;
  endo;
  retp(x);
endp;
