Last change
on this file since 37 was
37,
checked in by (none), 14 years ago
|
Added original make3d
|
-
Property svn:executable set to
*
|
File size:
1.5 KB
|
Rev | Line | |
---|
[37] | 1 | function y = trigamma(x) |
---|
| 2 | %TRIGAMMA Trigamma function. |
---|
| 3 | % TRIGAMMA(X) returns trigamma(x) = d**2 log(gamma(x)) / dx**2 |
---|
| 4 | % If X is a matrix, returns the trigamma function evaluated at each element. |
---|
| 5 | |
---|
| 6 | % Reference: |
---|
| 7 | % |
---|
| 8 | % B Schneider, |
---|
| 9 | % Trigamma Function, |
---|
| 10 | % Algorithm AS 121, |
---|
| 11 | % Applied Statistics, |
---|
| 12 | % Volume 27, Number 1, page 97-99, 1978. |
---|
| 13 | % |
---|
| 14 | % From http://www.psc.edu/~burkardt/src/dirichlet/dirichlet.f |
---|
| 15 | |
---|
| 16 | small = 1e-4; |
---|
| 17 | large = 8; |
---|
| 18 | c = pi^2/6; |
---|
| 19 | c1 = -2.404113806319188570799476; |
---|
| 20 | b2 = 1/6; |
---|
| 21 | b4 = -1/30; |
---|
| 22 | b6 = 1/42; |
---|
| 23 | b8 = -1/30; |
---|
| 24 | b10 = 5/66; |
---|
| 25 | |
---|
| 26 | % Initialize |
---|
| 27 | y = zeros(size(x)); |
---|
| 28 | |
---|
| 29 | % illegal values |
---|
| 30 | i = find(isnan(x) | (x == -inf)); |
---|
| 31 | if ~isempty(i) |
---|
| 32 | y(i) = nan; |
---|
| 33 | end |
---|
| 34 | |
---|
| 35 | % zero or negative integer |
---|
| 36 | i = find((x <= 0) & (floor(x)==x)); |
---|
| 37 | if ~isempty(i) |
---|
| 38 | y(i) = Inf; |
---|
| 39 | end |
---|
| 40 | |
---|
| 41 | % Negative non-integer |
---|
| 42 | i = find((x < 0) & (floor(x) ~= x)); |
---|
| 43 | if ~isempty(i) |
---|
| 44 | % Use the derivative of the digamma reflection formula: |
---|
| 45 | % -trigamma(-x) = trigamma(x+1) - (pi*csc(pi*x))^2 |
---|
| 46 | y(i) = -trigamma(-x(i)+1) + (pi*csc(-pi*x(i))).^2; |
---|
| 47 | end |
---|
| 48 | |
---|
| 49 | % Small value approximation |
---|
| 50 | i = find(x > 0 & x <= small); |
---|
| 51 | if ~isempty(i) |
---|
| 52 | y(i) = 1./(x(i).*x(i)) + c + c1*x(i); |
---|
| 53 | end |
---|
| 54 | |
---|
| 55 | % Reduce to trigamma(x+n) where ( X + N ) >= large. |
---|
| 56 | while(1) |
---|
| 57 | i = find(x > small & x < large); |
---|
| 58 | if isempty(i) |
---|
| 59 | break |
---|
| 60 | end |
---|
| 61 | y(i) = y(i) + 1./(x(i).*x(i)); |
---|
| 62 | x(i) = x(i) + 1; |
---|
| 63 | end |
---|
| 64 | |
---|
| 65 | % Apply asymptotic formula when X >= large |
---|
| 66 | i = find(x >= large); |
---|
| 67 | if ~isempty(i) |
---|
| 68 | z = 1./(x(i).*x(i)); |
---|
| 69 | y(i) = y(i) + 0.5*z + (1.0 + z.*(b2 + z.*(b4 + z.*(b6 + z.*(b8 + z.*b10))))) ./ x(i); |
---|
| 70 | end |
---|
Note: See
TracBrowser
for help on using the repository browser.