help-octave
[Top][All Lists]

## digamma & trigamma functions

 From: Christoph Mecklenbraeuker Subject: digamma & trigamma functions Date: Wed, 5 Jul 1995 20:18:01 +0200 (MET DST)

```Hi Guys,

A few days ago, I programmed the 'Digamma-' (aka 'Psi-') and 'Trigamma-'
functions which are part of the more general family of 'Polygamma-Functions'.
{which are recursively defined as derivatives of log(gamma(x))}

The code is tested  both for octave-1.1.0 and Matlab-4.2 on SunOS platform.

Caveat: Accuracy of the evaluation for *any* real argument was not my main
concern, because I need the functions only for integer and half-integer
arguments,
where *exact* finite recursions are readily availlable. So, I only
guarantee full precision for numbers of the form 1/2, 2/2, 3/2, 4/2,
5/2, 6/2, 7/2, ... but 3-4 decimals should be correct anywhere on the
real line. Maybe in near future, I will post an update with increased
accuracy.

So, here you get my alpha-version of both functions: suggestions for
improvements are most welcome. I tried to make the code as time
efficient as possible, but a lot of work still needs to be done.

Have fun!
Christoph

--(cut-here)--------------------file:digamma.m----------------------
function psi = digamma(z,method,debug)
%
%  psi = digamma(z)   ... Digamma-Function for real argument z.
%
%  digamma(z) = d/dz log(gamma(z)) = gamma'(z)/gamma(z)
%
%  if 'z' is a matrix, then the digamma-function is evaluated for
%  each element. Results may be inaccurate for real arguments < 10
%  which are neither integers nor half-integers.
%
%  psi = digamma(z,method)
%
%  possible values for optional argument 'method':
%    method = 1           : quick asymptotic series expansion (approximate)
%    method = 2           : finite recursion for integer values (exact)
%    method = 3           : finite recursion for half-integer values (exact)
%    method = 4 (default) : automatic selection of 1,2 or 3 for individual
%                           elements in z whichever is appropriate.
%

%  reference: Abramowitz & Stegun, "Handbook of Mathematical Functions"
%             Chapter "Gamma Function and Related Functions" :
%  implemented by: Christoph Mecklenbraeuker
%             (email: address@hidden), July 1, 1995.

dim = size(z);                          % save original matrix dimension
z = reshape(z,dim(1)*dim(2),1);         % make a column vector
I1 = ones(length(z),1);                 % auxiliary vector of ones

if(nargin==1)
method=4; debug=0;
elseif(nargin==2)
debug=0;
end;

if(debug == 1)                          % if debug==1: track recursion
[m,n] = size(z);
fprintf(1,'digamma: method = %d, size(z)=[%d %d],\t min(z)=%f,
max(z)=%f\n',...
method,m,n,min(min(z)),max(max(z)));
end;

if(method==1)                           % use 8th order asymptotic expansion
if(any(z<1))
fprintf(1,'Warning: some elements in argument of "digamma(z,1)" are < 1\n');
fprintf(1,'minimal argument = %g: digamma-result is
inaccurate!\n',min(min(z)));
end
% calculate powers of 1/z :
w1 = 1./z; w2 = w1.*w1; w4 = w2.*w2; w6 = w2.*w4; w8 = w4.*w4;
% generate coefficients of expansion: matrix with constant columns
a = [ -I1/2   -I1/12   I1/120  -I1/252  I1/240 ];
% make vector of powers of 1/z:
w = [  w1      w2       w4       w6       w8   ];
% calculate expansion by summing the ROWS of (a .* w) :
psi = log(z) + sum((a.*w).').';
elseif(method==2)
zmax = max(max(floor(z)));
psitab = zeros(zmax,1);
psitab(1) = -0.5772156649015328606;
for n=1:zmax-1;
psitab(n+1) = psitab(n) + 1/n;      % generate lookup table
end;
psi = psitab(z);
elseif(method==3)
zmax = max(max(floor(z)));
psitab = zeros(zmax+1,1);
psitab(1) = -0.5772156649015328606 - 2*log(2);  % = psi(1/2)
for n=1:zmax;
psitab(n+1) = psitab(n) + 2/(2*n-1); % generate lookup table
end;
psi = psitab(z+0.5);
elseif(method==4)                       % decide here which method to use
Less0 = find(z<0);                    % negative arguments evaluated by
reflexion formula
Less1 = find(z>0 & z<1);              % values between 0 and 1.
fraction = rem(z,1);                  % fractional part of arguments
f2 = rem(2*fraction,1);
Integers = find(fraction==0 & z>0);   % Index set of positive integer
arguments
NegInts  = find(fraction==0 & z<=0);  % Index set of positive integer
arguments
HalfInts = find(abs(fraction-0.5)<1e-7 & z>0); % Index set of positive
half-integers
Reals    = find(f2>1e-7 & z>1);       % Index set of all other arguments > 1
if(~isempty(Reals)) psi(Reals)    = digamma(z(Reals),1,debug);    end;
if(~isempty(Less1)) psi(Less1)    = digamma(z(Less1)+2,1,debug) - ...
1./z(Less1)-1./(z(Less1)+1);end;
% reflexion formula:
if(~isempty(Less0)) psi(Less0) = digamma(1-z(Less0),1,debug) -
pi./tan(pi*z(Less0)); end;
if(~isempty(Integers)) psi(Integers) = digamma(z(Integers),2,debug); end;
if(~isempty(HalfInts)) psi(HalfInts) = digamma(z(HalfInts),3,debug); end;
if(~isempty(NegInts))  psi(NegInts)  = Inf * NegInts; end;
end

psi = reshape(psi,dim(1),dim(2));

return;
--(cut-here)-----------------end-of-file:digamma.m----------------------

and here is the second one:
--(cut-here)--------------------file:trigamma.m----------------------
function y = trigamma(z,method,debug)

%  y = trigamma(z)   ... Trigamma-Function for real positive z
%
%  trigamma(z) = (d/dz)^2 log(gamma(z)) = d/dz digamma(z)
%
%  if 'z' is a matrix, then the digamma-function is evaluated for
%  each element. Results are inaccurate for real arguments < 10 which are
%  neither integers nor half-integers.
%
%  y = trigamma(z,method)
%
%  possible values for optional argument 'method':
%    method = 1           : quick asymptotic series expansion (approximate)
%    method = 2           : finite recursion for integer values (exact)
%    method = 3           : finite recursion for half-integer values (exact)
%    method = 4 (default) : automatic selection of 1,2 or 3 for individual
%                           elements in z whichever is appropriate.
%

%  reference: Abramowitz & Stegun, "Handbook of Mathematical Functions"
%             Chapter "Gamma Function and Related Functions" :
%  implemented by: Christoph Mecklenbraeuker
%             (email: address@hidden), July 4, 1995.

dim = size(z);                          % save original matrix dimension
z = reshape(z,dim(1)*dim(2),1);         % make a column vector
I1 = ones(length(z),1);                 % auxiliary vector of ones

if(nargin==1)
method=4; debug=0;
elseif(nargin==2)
debug=0;
end;

if(debug == 1)                          % if debug==1: track recursion
[m,n] =size(z);
fprintf(1,'trigamma: method = %d, size(z)=[%d %d],\t min(z)=%f,
max(z)=%f\n',...
method,m,n,min(min(z)),max(max(z)));
end;

if(method==1)                           % use 9th order asymptotic expansion
if(any(z<1))
fprintf(1,'Warning: some elements in argument of "trigamma(z,1)" are <
1\n');
fprintf(1,'minimal argument = %g: trigamma-result is
inaccurate!\n',min(min(z)));
end

% calculate powers of 1/z :
w1 = 1./z; w2 = w1.*w1; w3 = w1.*w2; w5 = w2.*w3; w7 = w2.*w5; w9 = w2.*w7;
% generate coefficients of expansion: matrix with constant columns
a = [ I1   I1/2   I1/6  -I1/30  I1/42 -I1/30];
% make vector of powers of 1/z:
w = [ w1   w2     w3     w5      w7    w9];
% calculate expansion by summing the ROWS of (a .* w) :
y = sum((a.*w).').';
elseif(method==2)
zmax = max(max(floor(z)));
ytab = zeros(zmax,1);
ytab(1) = pi^2/6;                     % = psi'(1)
for n=1:zmax-1;
ytab(n+1) = ytab(n) - 1/n^2;        % generate lookup table
end;
y = ytab(z);
elseif(method==3)
zmax = max(max(floor(z)));
ytab = zeros(zmax+1,1);
ytab(1) = pi^2/2;                     % = psi'(1/2)
for n=1:zmax;
ytab(n+1) = ytab(n) - 4/(2*n-1)^2; % generate lookup table
end;
y = ytab(z+0.5);
elseif(method==4)                       % decide here which method to use
Less0 = find(z<0);                    % negative arguments evaluated by
reflexion formula
Less1 = find(z>0 & z<1);              % values between 0 and 1.
fraction = rem(z,1);                  % fractional part of arguments
f2 = rem(2*fraction,1);
Integers = find(fraction==0 & z>0);   % Index set of positive integer
arguments
NegInts  = find(fraction==0 & z<=0);  % Index set of positive integer
arguments
HalfInts = find(abs(fraction-0.5)<1e-7 & z>0); % Index set of positive
half-integers
Reals    = find(f2>1e-7 & z>1);       % Index set of all other arguments > 1
if(~isempty(Reals)) y(Reals)    = trigamma(z(Reals),1,debug);    end;
if(~isempty(Less1)) y(Less1)    = trigamma(z(Less1)+2,1,debug) + ...
1./z(Less1).^2+1./(z(Less1)+1).^2;end;
% reflexion formula:
if(~isempty(Less0)) y(Less0)=
-trigamma(1-z(Less0),1,debug)+(pi./sin(pi*z(Less0))).^2; end;
% integers:
if(~isempty(Integers)) y(Integers) = trigamma(z(Integers),2,debug); end;
% half-integers:
if(~isempty(HalfInts)) y(HalfInts) = trigamma(z(HalfInts),3,debug); end;
% negative integers:
if(~isempty(NegInts))  y(NegInts)  = Inf * NegInts; end;
end

y = reshape(y,dim(1),dim(2));
return;
--(cut-here)--------------------end-of-file:trigamma.m----------------------

--
Christoph Mecklenbraeuker             | send electronic mail to:
Lehrstuhl f. Signaltheorie (IC 5/35)  | address@hidden
Ruhr Universitaet Bochum              | Tel: +49(234) 700 6119
D-44780 Bochum, Germany               | Fax: +49(234) 709 4261
---
```